Grizzled Skippers stuck in the south: Population‐level responses of an early‐successional specialist butterfly to climate across its UK range over 40 years

Climate change has been predicted to facilitate poleward expansion of many early‐successional specialist invertebrates. The Grizzled Skipper, Pyrgus malvae, is a threatened butterfly in long‐term decline that has not met expectations of northern expansion in Britain, possibly indicating that climate change has not improved northern habitat suitability or that another driver (e.g. land use change) is masking its effects. Here, we explore the effect of climate on population size trends over four decades, and whether any regions show an improving population trend that may be a precursor to northern expansion. Examining detailed spatio‐temporal abundance data can reveal unexpected limitations to population growth that would not be detectable in widely used climate envelope models.


| INTRODUC TI ON
Temperatures have been increasing globally since the mid-20th century (Intergovernmental Panel on Climate Change, 2014; Met Office, 2019), with increasing regularity of extreme heat and drought summer events, and longer growing seasons (Met Menzel et al., 2020;Office, 2019). It has been predicted that in a scenario of ambient warming, species' geographical ranges may shift in order to track their climate envelopes, predominantly polewards and to higher altitudes (Parmesan et al., 1999;Root et al., 2003). For earlysuccessional specialists that frequently require warmer microclimates within their range, range changes may result from an ability to occupy later seral stages under warming conditions, hereby becoming less specialized and less reliant on early-successional habitats (Komonen et al., 2004;O'Connor et al., 2014;Thomas et al., 1999).
Less specialized habitat requirements may enable species to expand their cover within and beyond current population boundaries (Settele et al., 2008;Thomas et al., 2001). Some early-successional species, such as the Brown Argus butterfly (Aricia agestis), appear to have already broadened their habitat niche as a result of warmer summers, and expanded polewards (Pateman et al., 2012). However, one early-successional specialist that has not responded to climate change as expected is Pyrgus malvae, the Grizzled Skipper butterfly, that has declined in the UK by 55% between 1976 and 2018 (Brereton et al., 2019).
Expectations of range expansions and shifts are not unique to early-successional specialists (Chen et al., 2011;Davis & Shaw, 2001), but in particular P. malvae, like other southern restricted and warmpreferring species, is expected to benefit from warming temperatures in large part due to its requirements for warm microclimates for oviposition and development during temperature-sensitive immature stages (Krämer et al., 2012;Weiss et al., 1988). As earlysuccessional specialists like P. malvae often rely on the presence of warm microclimates to persist in cooler areas, such as at higher altitudes or more poleward locations, their habitat breadth in those cooler locations may be narrower than in warmer areas in southern Europe . Therefore, an increase in temperature in the northern parts of P. malvae's range may increase breeding habitat availability and facilitate increases in abundance and subsequent range expansions. However, contrary to predictions of poleward expansion under a warming climate, no indication of natural northern expansion has been observed in P. malvae (National Biodiversity Network 2020), and declines have been observed that on average do not appear to have slowed in the last decade, despite warming conditions (Brereton et al., 2019;Fox et al., 2015). Similar observations of species not fully tracking their expanding range limits have been made across several taxonomic groups and countries (Bertrand et al., 2016;Devictor et al., 2012;Vanderwal et al., 2013;Warren et al., 2001).
Although there has not yet been evidence of northern expansion by way of colonizations beyond the current north range edge, investigating changes in abundances of northern populations could provide insight into the likelihood of imminent expansion. Increasing abundances close to the range edge of Pyrgus malvae have been found previously to represent an intermediate step prior to colonization (Maggini et al., 2011). Therefore, investigating whether P. malvae abundance is showing increases at poleward range limits over time could provide an early indication of increasing climate suitability as would be expected from spatial-based predictions. Spatial-based predictions alone may not be able to fully capture limitations to expansion, such as dispersal capacity and biotic habitat requirements, so assessing trends over time and space may offer a more accurate reflection of potential range changes. In addition, exploring the effects of climate variables on population abundance could provide further insight into the role that climate plays in influencing longterm trends, and differences in regional trends if they are found.

| Aims
Here, we aim to explore the effect of climate on long-term population size trends, and whether UK Butterfly Monitoring Scheme (UKBMS) indices of population size at UKBMS sites are decreasing at differing rates between regions. Population size trends differing with latitude/longitude is of interest as relatively positive trends in population size may be a precursor to northern expansion (Maggini et al., 2011;Mair et al., 2014).
We hypothesize that if climate change is ameliorating declines in P. malvae in the UK, climate predictors will account for a large portion of variation in abundance otherwise explained by time, thereby emerging as a dominant driver of long-term trends. In addition, if climate is the limiting factor for northern range expansion, but range edge conditions are becoming more suitable under climate change, abundance will have relatively positive long-term trends at higher latitudes.

| ME THODS
Pyrgus malvae is found in early seral stages of grassland, woodland and brownfield habitats (Brereton et al., 1998) in the UK. However, its range has contracted, having lost its northern most population in York in the 1990s and having gone locally extinct in several counties, including in Derbyshire in 2007, where it is currently the focus of reintroduction efforts (Butterfly Conservation, 2018). Another reintroduction was attempted at Gait Barrows in Lancashire in 2002 (geographical locations shown for reference in Figure 1), where it had occupied habitat in the early 20th century, but the population is believed to have gone extinct in 2007 (Coleman, 2017).
The UK Butterfly Monitoring Scheme is a long-term coordinated recording scheme for butterflies in which a core component is a network of more than 2000 sites at which adult butterflies are counted weekly over the flight season between April and September (Brereton et al., 2019;Pollard & Yates, 1993). The weekly butterfly counts are used to produce annual indices that account for missing weeks, and this index can be used as a proxy for a species' population size at a site . Henceforth, "Population Index" is used to mean the P. malvae annual abundance index at a given UKBMS transect site, which is calculated from a generalized additive model . Mixed models were used to in- Climate data were obtained from Met Office data at 5 × 5 km resolution, as monthly mean temperature, and monthly rainfall values (mm) (Perry, 2005). Each population index value was matched to the climate data from its own 5 km square, for the months from June of the previous year to May of the year of the record (inclusive).
This time window was chosen to span the life stages of a single generation from the larval period to the following flight period. Mean monthly temperature and rainfall were chosen for ease of interpretation and generality (Palmer et al., 2015;Pollard, 1988). Preliminary investigation of three-month pooled climate values did not produce notably different or more significant effects on population index than monthly mean temperature and rainfall (Appendix S1).
A total of 483 sites, mapped in Figure 1, were available that had at least one year with sufficient data to produce an annual population index during the survey period, and fell within a 5 km grid square with Met Office data. Year-site combinations with P. malvae presence but with too few visits to estimate a population index (according to UKBMS criteria) were omitted.
Generalized linear mixed effects models were run using the R package "lme4" (Bates et al., 2015). The error family used was negative binomial to account for a high proportion of zeros in the population index data (38%). To prevent convergence issues, all predictor variables were standardized by subtracting the mean of each variable, and dividing by the standard deviation.
In all models, to account for pseudoreplication, as sites were visited repeatedly, and for stochastic variation between years and regions across the country, year, site number and 50 × 50 km grid squares (Total: 59) were introduced as random factors. Our approach did not account for population size in the previous year, because of the limited number of UKBMS records with visits to the same site the prior year. 50 × 50 km square was selected as the appropriate size in order to minimize the number of grid squares with very few sites, and at a high enough resolution to isolate geographic aggregations of sites with unusually high residual errors from the model when 50 × 50 km square is not included as a random effect. Model fit was evaluated using AIC comparisons, and validated by assessing heterogeneity in residuals.
We ran three sets of statistical models to elucidate the effects of climate versus potentially co-varying effects of geography. The final minimum adequate model was then used to explore the impact of climate on long-term population index trends.

| Non-climate predictors
To assess how population size has changed over time across the UK geographical range, irrespective of climate, models were first run containing only year, latitude and longitude as fixed effects. A model was run for each of these fixed effects independently, and then with F I G U R E 1 Map showing the geographical distribution of UKBMS sites included within analyses, where each site is represented by a black circle. The black dashed lines depict the total range of all UKBMS sites with Grizzled Skipper records. The red dashed lines depict the median latitude and longitude of all records used in analyses. Locations of P. malvae extinctions and reintroductions are shown as red symbols; the Derbyshire reintroduction site (cross), Lancashire reintroduction site (triangle), Yorkshire extinction site (diamond) all effects, and two-way interactions, with AIC stepwise regression to obtain the minimum adequate model (MAM01). The initial maximal model (MM01) therefore can be described as: Upon running the final model, we observed a large difference between the real and predicted population index values at the beginning of the recording period in the north-west of the study area.
As we believed this to be due to data scarcity for this space and time, we reran a model omitting data prior to 1990, and the outputs were qualitatively unchanged (Appendix S2).

| Climate predictors
To explore the effect of climate in each month independently, and to avoid auto-correlation between months, individual models were run for each monthly climate value as the only fixed effect, for both temperature and rainfall (24 total models).
An example of one such model can be described as: All variables found to have a significant effect were included in a single model, with no interactions, to better identify variables that may be accounting for similar variation in population index. The least significant fixed effect was removed in a stepwise manner to minimize AIC and arrive at the minimum adequate model for climate effects (MAM02).

| All predictors
To assess how much of the trends over time and geography ob- The initial maximal model (MM03) can be described as:

| Predicting indices under different climate conditions
Having established the minimum adequate model containing both climate and non-climate predictors (MAM03), we explored the extent to which fixing local (5 × 5 km) climate conditions at the mean climate values for the first five years of the study (1976)(1977)(1978)(1979)(1980)  In addition, to further investigate possible trends in climate suitability for P. malvae over space and time, we estimated climate suitability using just climate effect coefficients from MAM03 and mean climate values for four regions within the study area, the dimensions of which are presented in Figure 1. The outputs are presented in Appendix S4.

| RE SULTS
Population index henceforth refers to the P. malvae annual abundance index at a given UKBMS transect site. For models involving only non-climate predictors (Methods Section i.), year, in isolation, had a significant negative effect on population index (Est. coef = −0.16, Std. error = 0.05, p = .001), but no effect was observed of latitude (Est. coef = −0.19, Std. error = 0.11, p = .088) or longi- Population index ∼ Previous July Temperature  (Table 1, geographical differences illustrated in Figure 3).
When considering effects of climate one by one (Methods Section ii.), we found significant negative effects of December and January temperature, and previous June and previous July rainfall, and a positive effect of February rainfall (Appendix S1). Not all of these effects remained significant when added to a single model; the effects of January temperature, February rainfall and previous July rainfall were the least robust (Appendix S1; Table 1).
The fit of the minimal adequate non-climate model (MAM01) was improved by the inclusion of December temperature, and February, June, and July rainfall effects (ΔAIC from MAM01 = −20.2, Table 1).
Interestingly, neither July nor February rainfall effect was significant when both were included in the model, but when either variable was singly removed, the remaining variable became significant, and the fit of the model was slightly improved (Table 1) The inclusion of climate slightly reduced the effect size of longitude and latitude on trends over time (Table 1), but the effect size of climate itself remained small overall (Figure 2). In the final model (MAM03), a one degree increase in December temperature was estimated to account for a log decrease of 0.08 in population index (Table 1). Similarly, a 5cm increase in June rainfall accounted for a log decrease of 0.09 in population index. In contrast, a 5cm increase in February rainfall accounted for a log increase of 0.08 in population index ( Table 1).
The effect size of climate is shown in Figure 2

| D ISCUSS I ON
In this investigation, we explored spatial, temporal and climatic effects at a 5-km resolution on P. malvae population size at monitored UKBMS sites, in order to increase understanding of the role of geography and climate in the long-term trends of the declining butterfly. We used a mixed-model approach to control for increasing survey effort over time, and stochastic variation between years and regions. We observed no significant positive effects of temperature on population index; population index henceforth refers to the P. malvae annual abundance index at a given UKBMS transect site. In the best fitting model, December temperature and June rainfall had significant negative effects on population index, while February rainfall had a positive effect on population index. However, the effect sizes of all climate variables were very small TA B L E 1 The coefficients of fixed effects, and model weightings from a subset of models run with the same response data

| Climate effects
We found no significant positive impacts of temperature on population index in any month (Appendix S1), which would not suggest that the P. malvae is overall benefitting from climate warming during its temperature-sensitive life stages.
To the contrary, we observed a negative effect of December temperatures on population index, which corroborated established detrimental effects of warm winters on a variety of invertebrates (McDermott Long et al., 2017;WallisDeVries & van Swaay, 2006).
Suggested mechanisms for which include higher fungal infection rates (Harvell et al., 2002;Radchuk et al., 2013), and phenological asynchrony (Parmesan, 2007)  However, early-successional specialists, such as P. malvae, frequently rely on warm microclimates even when immature life stages occur after spring (Brereton et al., 1998;de Schaetzen et al., 2018;Asher et al., 2001), particularly in cooler locations (e.g. at cool range margins). Early-successional specialists may therefore also be impacted by greater total vegetative biomass, and cooler microclimates, resulting from longer and earlier growing seasons following warmer winters. Similarly, populations of temperature-sensitive species in general at their cool range edges could be affected by cooling of microclimates relied upon for refuge from cooler ambient conditions (Oliver et al., 2009).
Lusher vegetation growth may also be a potential mechanism underlying the observed negative effects of June rainfall on P. malvae abundance (Appendix S1), because late instar caterpillars are feeding in June; similar effects have been observed previously in P.
Given that a model not including the effect of February rainfall was of a similar fit to the minimal adequate model with February rainfall (MAM03), we are cautious about the implications of this effect. The reported effect of February rainfall may be an artefact of its correlation with July rainfall in this data set, given a negative effect of summer rainfall on P. malvae abundance having precedence in the literature (Pollard, 1988;Roy et al., 2001;WallisDeVries et al., 2011).

| Role of climate in declines
While the climate effects were significant, the effect sizes were small, and they did not account for long-term trends in population index, or for geographical differences in those trends. Modelling population index over time revealed similar long-term trends in both unchanging, and observed climate scenarios (Figure 2). Climate's limited role in long-term trends may result from an absence of significant trends in combined climate suitability between 1976 and 2016 in P. malvae's UK range (Appendix S4). Climate did appear to play a larger role in short term yearly variation (Figure 2). Similar findings regarding climate driving short-term variation in abundance between years in a related species, Pyrgus armoricanus, were also recently reported, in conjunction with a limited effect of temperature on colonization and extinction dynamics relative to the effects of connectivity and host plant density (Fourcade et al., 2017).
Furthermore, inclusion of climate predictors in the final model did not qualitatively change the effects of interactions between geography and time on P. malvae population index (Table 1), suggesting the climate effects are not the drivers of steeper temporal declines in the north and west than the south and east, respectively. Correspondingly, climate suitability does not show differing temporal trends in any of the geographic subsets shown in Figure 3 (Appendix S4).
Although this study's scope does not extend to estimating P.

| Drivers of decline
We show strong evidence of long-term abundance declines between 1976 and 2016, especially in the northern and western parts of the UK range. We argue that any simple climate mechanism behind these declines would have been identified with our statistical tests and therefore are left with the conclusion that a non-climatic driver is behind the declines, most likely habitat loss and degradation (Brereton et al., 1998).
Furthermore, the quality of remaining habitat has deteriorated in places due to abandonment and changes in natural grazing pressure, such as those resulting from losses of rabbit populations to myxomatosis (Sumption & Flowerdew, 1985;Travers et al., 2019).
In recent decades, some semi-natural grasslands have been restored or improved through rabbit population recoveries , and agri-environment scheme introductions, with some proven benefits to butterflies (Brereton et al., 2011;Davies et al., 2005). However, early-successional deciduous woodland habitat quality continues to show an overall decline in the UK, owing to management costs and logistic difficulties. (Atkinson & Townsend, 2011;Keith et al., 2009;Kirby et al., 2017). Contrastingly, coppiced woodland status varies between regions of Central Europe, as a result of silvicultural intensification and abandonment (Bergmeier et al., 2010), leaving the future status of P. malvae's habitat availability uncertain.
Additionally, north-west populations may have been more sensitive to changes in habitat quality than those in the southeast, particularly at the historical northern cool range margin.
Populations at cool range margins can be more reliant on warm microclimates than those in the centre (Hodgson et al., 2015;Suggitt et al., 2018), so shifts to later successional stages may be more detrimental to cool range edge-based populations. P. malvae has also been shown to have higher habitat specificity at sites with higher rainfall, potentially amplifying detrimental effects of succession (Mair et al., 2012).
This brings us back to the idea that climate warming (in general, but especially in winter and spring) can lead to accelerated vegetation growth in temperate climates and reduce bare ground availability, essentially accelerating succession (Zhou et al., 2001). Some early-successional species may be affected more by this than by the opposing effect that ambient warming makes even tall vegetation warmer. This makes it particularly difficult to plan for the conservation of early-successional habitat specialists, which tend to be threatened across Europe (Kuussaari et al., 2007;van Swaay et al., 2006) as well as in the UK (Thomas, 1993;Thomas et al., 2015). We might expect each species to respond idiosyncratically to climate change at its cool range margin, depending on the precise balance of the opposing effects.
Habitat availability and climate suitability interact, so for species like P. malvae, it is a conservation priority to boost populations by improving habitat condition, both for the direct benefits obtained, and to indirectly allow the species to track climate if and when suitability improves. There is increasing empirical evidence that climatedriven range shifts are slower or impossible if there is insufficient habitat availability (Devictor et al., 2012;Platts et al., 2019). For P.
malvae, overall declines in population size over time, and more negative trends in the north than the south, as was observed (Figure 3), suggest a low likelihood of poleward expansion in the near future (Maggini et al., 2011). Nevertheless, alleviating the non-climatic threats to the species could lead to an expansion of range as well as abundance. This could be important for the global status of the species because as the entirety of the UK is still northerly within P.
malvae's overall range, and it may be suffering from heatwaves or droughts at its extreme southern limits.
The complex interactions between habitat and climate discussed so far complicate predictions of habitat suitability when considering range shifts. Species distribution models based only on spatial data (e.g. occurrence) may not fully account for other habitat requirements that limit the potential for a species to exploit opportunities in otherwise climatically suitable habitats. Predictions of range shifts/ expansions are particularly vulnerable in a non-equilibrium system, such as one undergoing changes in climate and land management as is the case in many systems across the globe (Elith & Leathwick, 2009).

DATA ACCE SS I B I LIT Y S TATE M E NT
The data that support the findings of this study are available upon

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13245.

B I OS K E TCH
Fiona Bell's current research focus is on climate change, land use change, the interactions of drivers of decline, and the conservation of threatened species, with a particular interest in lepidoptera, that is the Grizzled Skipper (Pyrgus malvae), which is the focal species of her ongoing PhD.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.