Summer drought decreases the predictability of local extinctions in a butterfly metapopulation

The ecological impacts of extreme climatic events on population dynamics and/or community composition are profound and predominantly negative. Here, using extensive data of an ecological model system, we test whether predictions from ecological models remain robust when environmental conditions are outside the bounds of observation. First, we report a 10-fold demographic decline of the Glanville fritillary butterfly metapopulation on the Åland islands (Finland). Next, using climatic and satellite data we show that the summer of 2018 was an anomaly in terms of water balance and vegetation productivity indices across the habitats of the butterfly, and demonstrate that population growth rates are strongly associated with spatio-temporal variation in climatic water balance. Finally, we demonstrate that covariates that have previously been identified to impact the extinction probability of local populations in this system are less informative when populations are exposed to (severe) drought during the summer months. Our results highlight the unpredictable responses of natural populations to extreme climatic events.


INTRODUCTION
One of the major challenges in conservation biology is to identify the species and populations that are most vulnerable to extinction. Long-term monitoring of ecological model systems, both at local and global scales, have facilitated conservation objectives by identifying the factors affecting population declines and extinctions (Hanski et al. 1995;Pimm et al. 2014;Dornelas et al. 2019). In addition, ecological models aiming to improve our understanding of population dynamics in temporally varying environments have been employed to shed light on which regions should receive priority for conservation and to predict which species and/or populations are most vulnerable (e.g. Franklin et al. 2014;Oliver et al. 2015). In temporally autocorrelated environments, where conditions tend to change gradually, these predictions may be reliable over short timescales, and can therefore be used to make conservation efforts more effective.
When populations are exposed to conditions that are beyond the normal range, such as in the case of extreme climatic events (ECEs), the factors underlying population dynamics may be less relevant and consequently predictions relying on these factors less reliable. ECEs have increased in recent decades, not only in frequency but also in intensity and duration (Jentsch et al. 2007;Coumou & Rahmstorf 2012;Ummenhofer & Meehl 2017). Recent studies have demonstrated that even a single ECE, such as a flood or a drought, can have profound impacts on natural populations (Bailey & van de Pol 2016;Altwegg et al. 2017;Grant et al. 2017), and can even cause the collapse of an entire ecosystem Harris et al. 2018). For example, a pan-tropical episode of coral bleaching, triggered by a marine heatwave in 2016, eradicated over 60% of the coral communities in the Great Barrier Reef (Hughes et al. 2017).
Studies assessing impacts of ECEs are generally conducted by performing perturbation experiments (e.g. Bokhorst et al. 2008;Pansch et al. 2018) or by opportunistically taking advantage of a rare event (e.g. Smith 49 50 51 52 (ongoing) systematic monitoring, the system provides a unique opportunity to improve our understanding of how large spatially-structured natural populations respond to extreme climatic conditions, and to study how ECEs affect the value of predictive models.
Since 1993, the occupancy and abundance of local populations of the butterflies, across a network of about 4400 dry meadows, has been systematically monitored. Estimates of local population sizes are acquired during autumn when all potential habitat patches are surveyed for the presence of the winter nests (for details see; Ojanen et al. 2013). Despite increasing year-to-year fluctuations (Hanski & Meyke 2005;Tack et al. 2015), the overall size of the metapopulation has remained relatively stable, with about 20% of the available habitat patches being occupied each year. Survey data have been used to demonstrate that the long-term persistence and population dynamics (local extinctions and re-colonizations) of the metapopulation largely depends on the number of available habitat patches, their area, and their spatial configuration (i.e. connectivity; Hanski & Ovaskainen 2000;Dallas et al. 2019). Other key factors impacting local extinctions and population dynamics are the habitat quality (e.g. amount of host plants; Harrison et al. 2011;Schulz et al. 2019) and genetic characteristics of the individuals (Saccheri et al. 1998;Niitepõld & Saastamoinen 2017).
Annual changes in population size across the metapopulation are mainly driven by the growth rates of local populations over summer (r = 0.96), rather than being associated to variation in overwinter mortality (r < 0.01; see Tack et al. 2015 for details). Moreover, we recently demonstrated that variation in population growth rates is strongly associated with variation in temperature and precipitation across the habitat patch network (Kahilainen et al. 2018). Spatial variability in these climatic conditions contributes to the overall stability of the metapopulation by ensuring that local extinctions are compensated by stable or increased population sizes and colonization rates in other areas of the metapopulation. As the climatic conditions have become more homogeneous across Åland in the last decade, the Glanville fritillary metapopulation has now become more spatially synchronous in its demographic dynamics (Kahilainen et al. 2018), potentially making it more vulnerable to large-scale environmental perturbations (Hanski & Woiwod 1993).
Here, we identify the key anomalies that occurred in the summer of 2018 on the Åland islands and demonstrate that exposure to these extremes was sufficient to drive the regional population declines of the butterfly observed across the metapopulation. Secondly, we explore how the ECE impacts the performance of ecological models by testing whether the covariates that have previously been identified to affect local extinctions in this system remain informative under extreme climatic conditions.

Long-term data survey
On the Åland islands, the Glanville fritillary inhabits dry meadows and pastures with at least one of the two larval host plant species, the ribwort plantain (Plantago lanceolata) or the spiked speedwell (Veronica spicata) present. The entire study region (50 × 70 km) has been mapped for potential habitat patches, with a total of 4,415 patches (September 2015). The population occupancy (presence of a larval nest within a habitat patch) and abundance (total number of larval nests within a habitat patch) is assessed annually when the larvae have entered the overwintering stage within their conspicuous larval nests. All the habitat patches on the archipelago are inspected with the help of approximately 50 field assistants. In the field, the locations of the larval nests are marked and data are stored into the EarthCape database (Ojanen et al. 2013). The occupancy and abundance data comprise raw counts, which is a function of both population size but also of detection probability. Estimates from control surveys conducted before 2018 suggest that the presence of the butterfly may be missed in up to 15% of the occupied patches, with non-detection mainly occurring in very small populations (Hanski et al. 2017). Both NDVI and EVI data are available on a 16-day basis for a 19-year period between 2000 and 2018. Data are derived from bands 1 and 2 of the MODIS on board NASA's Terra satellite. A time-series of NDVI and EVI observations can be used to examine the dynamics of the growing season, ecosystem health or monitor phenomena such as droughts. EVI is an 'optimized' vegetation index designed to enhance the vegetation signal with improved sensitivity in high biomass regions and improved vegetation monitoring through a decoupling of the canopy background signal and a reduction in atmosphere influences. NDVI and EVI values were calculated for each of the 4,415 habitat patches from April to August. Monthly temperature and precipitation values and water balance (precipitation -potential evaporation) are derived from Jomala precipitation minus potential evapotranspiration (PET). PET is the amount of evaporation and transpiration that would occur if a sufficient water source were available. PET was calculated using the Penman-Monteith equation (Guo et al. 2016). Based on the climatic water balance values a time-series of the Standardized Precipitation-Evapotranspiration Index (SPEI) was calculated (Vicente-Serrano et al. 2009). An important advantage of the SPEI is that it can be computed at different time scales ( i.e. it is possible to incorporate the influence of a variable range of past values in the computation). For example, a SPEI-3 value implies that data from the current month and those of the past two months were used to compute the drought index.

Population growth rate model
Following an analysis pipeline developed for the Glanville fritillary metapopulation (Kahilainen et al. 2018), we fitted a Bayesian linear mixed effects model to study the association between weather conditions and population growth rate at the regional scale. Instead of using temperature and precipitation values, we used monthly water balance indices, which are more explicitly associated with drought conditions. These values were computed for the whole study area at the spatial resolution of 1 km 2 using national climate station observations covering the period of 1993-2018. Additionally, national observations were supplemented by records from the surrounding countries, and kriging interpolation was applied to account for topography and water bodies (Aalto et al. 2016). We focused on the dynamics of groups of spatially clustered semiindependent habitat patch networks (SIN; see Hanski et al. 2017) instead of local populations, because populations in individual patches frequently go extinct, and estimates of population growth at this spatial scale would therefore be very heterogeneous. Furthermore, we chose to focus on a subset of 33 SINs that are considered viable in terms of the spatially explicit metapopulation theory and were occupied most of the study period (Hanski et al. 2017). Yearly population growth rates (r) were derived as: , where N SIN,t is the number of overwintering larval nest found in a SIN in the fall of year t. A single nest was added to both the numerator and the denominator to account for the occasional event that all patches in a SIN were unoccupied. We then extracted water balance data for each SIN (see previous section) and fitted a linear mixed model for population growth rate using the water balance values across different life-history stages of the butterfly as covariates (see also Radchuk et al. 2013). To reduce the number of variables in the model, but keeping those from months known to be of importance for population growth (see Tack et al. 2015;Kahilainen et al. 2018), an averaged water balance value (September-February) was included to reflect the conditions during the diapause period. In the model we included a random intercept and a first order autocorrelation term for each SIN. We implemented the model in Stan statistical modeling platform (Carpenter et al. 2017), using R packages brms (version 2.7.0; Bürkner 2017) and RStan (2.17.3; Stan Development Team 2018) as its interface. For further details regarding the analysis pipeline and implementation of the model, see Kahilainen et al. 2018. From the model we then obtained fitted annual growth rates for each of the 33 SINs included in the data and compared these with the observed growth rate values (Fig. 3).

Extinction probability model
We fitted a linear mixed effect model (logit link function) to binary patch extinction data calculated for each year of the survey between 1999 and 2018 to examine whether the high extinction rates in 2018 could be attributable to covariates previously recognized to impact the extinction probability of local populations (Hanski et al. 1995;Saccheri et al. 1998). Year was included as a random effect. Environmental variables included as fixed effects were the natural logarithm of patch area, the amount of the dominant host plant in each patch estimated on a scale of 1-3, the percent of host plants that were dried out, and the percent of the patch that was grazed (see Ojanen et al. 2013 for more details). Historical contingency was accounted for by considering the natural logarithm of the number of larval nests in each focal patch in the previous year, and the number of larval nests in the neighborhood of each focal patch in the previous year, defined as: where, d is the distance between focal patch i and j, and Nj is the number of nests in patch j in the previous year. We trained the model using 80% of the data, and subsequently assessed model performancequantified using the area under the receiver operating characteristic (AUC) -on the remaining 20% of the data. This procedure was repeated 10 times in order to assess the sensitivity of model performance to the train/test split.

Demographic declines across the butterfly metapopulation in 2018
Despite systematic surveying efforts, an all-time low of only 91 larval nests were recorded during the autumn survey in 2018 (the average number of nests recorded each year is 2750; figure 1). The number of occupied habitat patches was also an order of magnitude lower than in any average year (N=54; figure 1A), and only a single colonization event was recorded. The number of larval nests recorded within the occupied habitat patches was approximately 65% below the average (table S1).

Drought drives regional population declines
We examined whether the extreme climatic conditions experienced by the butterfly in 2018 were sufficient to drive population declines at the scale of the semi-independent patch networks (SINs; Hanski et al. 2017).
Overall, monthly water balance indices -especially the water deficits during May and July -were highly associated with population declines at the SIN level ( figure 3; table S2). With few exceptions, the model residuals for 2018 were negative, indicating that the observed declines in most SINs were more severe than suggested by the model fit ( figure 3). Indeed, SIN level population declines in 2018 were among the most dramatic ones observed since the start of the annual survey, with nearly 80% of the growth rate values falling within the lowest 5 percentile (figure 3; table S1). We first confirmed the previously identified factors affecting extinction probability of local populations in this system. In short, the probability of extinction is lower in habitat patches that are large and wellconnected. Patches with low population densities in the previous year, reduced host plant availability, or those used for livestock grazing, have higher extinction probabilities (table S3). Secondly, we estimated the predictive accuracy of our environmental extinction model. The mean AUC score across the 20 years was 0.76, implying that model covariates predict the extinction of a local population reasonably well (figure 4a; AUC scores above 0.80 reflect good discrimination (Swets 1988). The AUC scores varied significantly across years with the lowest values observed in years that are characterized by water deficits during summer (figure 4b). Model covariates were more likely to predict the extinction probability correctly in the years prior to the ECE of 2018.

DISCUSSION
The ongoing climate crisis is expected to result in a rapid increase in the frequency, intensity and duration of extreme climatic events (e.g. Sun et al. 2014;Fischer & Knutti 2015). Owing to the profound impacts these extreme events can have on ecological functioning, the vulnerability of natural populations to the impacts of climate change may be severely underestimated Ummenhofer & Meehl 2017).
The Åland islands, like many other regions across the northern hemisphere, sustained severe drought in the summer of 2018. Using extensive long-term data of the spatially structured natural population of the Glanville fritillary butterfly, we demonstrate that both the climatic conditions during the summer and biological response of the metapopulation exceeded the 5% threshold values typically used to define ECEs (van de Pol et al. 2017). The majority of studies investigating the impacts of ECEs have focused on the ability of organisms to cope with the direct effects of the extreme environmental conditions. For example, the overwintering life stages of some species of butterfly are known to be very sensitive to warm winter temperatures (Radchuk et al. 2013;Roland & Matter 2016). While an organisms' critical physiological limits are likely to be important risk factors (Hoffmann & Sgrò 2018), extinctions of local populations could also be triggered indirectly (Maron et al. 2015). Our estimates of the occupancy and abundance of local populations across the archipelago are collected in the autumn, and hence only represent populations that survived the potentially harsh summer conditions. Therefore, we cannot infer whether the population declines and local extinctions were triggered directly by the arid conditions, or indirectly through for example climate-induced restrictions in resource availability (Maron et al. 2015;Johansson et al. 2020).
Nonetheless, using satellite-derived vegetation indices we demonstrate here that productivity was substantially reduced within the habitat patches of the butterflies in 2018. These data suggest that the severe drought negatively impacted the quality and availability of the host plants of these herbivorous insects, making resources scarce when demands were high. These results are in line with detailed observations made in 12 focal patches in 2018, where despite observing large numbers of clutches in early summer only two larval nests were found in autumn (Salgado et al. under review). The host plants are known to recover quickly from drought, and the proportion of desiccated host plants in the habitat, as assessed during the survey in autumn, has been shown to be a poor predictor of plant condition during larval development (Tack et al. 2015).
We explored whether the strong deviations from usually experienced conditions influenced our ability to forecast the short-term biological response of the metapopulation, both at a regional and local scale. First, we demonstrate that the population growth rates at regional scales (i.e. the level of SINs; Hanski et al. 2017) are positively associated with climatic water balance during summer months (table S2). In other words, population sizes of this drought-sensitive butterfly species declined in years that were characterized by water deficits during summer (see also Tack et al. 2015;Kahilainen et al. 2018). This result strongly points towards climatic conditions as the key driver of the dramatic decline of the metapopulation in 2018. In addition, we find that the relationship between the climatic variables and the biological response of the metapopulation is largely linear (i.e. model residuals of 2018 are negative, but generally do not exceed model confidence limits; figure 3; van de Pol et al. 2017). These results highlight that our detailed understanding of the population dynamics of this ecological system allowed us to -at least in part -forecast the observed population declines in response to the described extreme event.
While we were able to capture climate-induced changes in dynamics at regional scales, we were less successful in predicting the extinction probabilities of local populations in 2018 (figure 4a). In general, the model covariates, such as habitat size and connectivity, were likely to predict the extinction probability correctly. In 2018, even populations inhabiting large and well-connected patches with a high abundance in the previous year went extinct. This suggests that summer droughts, and in particular the extreme climatic other than those highlighted in previous studies with the system (reviewed in Ovaskainen & Saastamoinen 2018). Hence, uncharacterized ecological variables, such as the water holding capacity of the patch and/or the microhabitat heterogeneity within each meadow or pasture, may potentially be important determinants of population persistence under arid environmental conditions. In addition, interactions with other species in the community (e.g. parasitoids or predators) may be different in dry years. Finally, since annual overwinter survival varies between 0.41 and 0.83 (Tack et al. 2015), extending the current models to include variables that capture variation in winter conditions could potentially improve their predictive value.
As a potential caveat we note that, despite using a systematic and methodical protocol, our probability to detect the presence of the butterfly is known to be reduced in very small populations (Hanski et al. 2017).
Hence, strong climate-induced declines in population sizes may therefore contribute to decreased detection rates (and, consequently, an overestimation of the biological response). Other variables, such as reduced nest sizes or shifts in habitat choice, may potentially also negatively affect detectability under climatic extremes.
However, our surveying protocol dictates that when a patch is judged to be unoccupied, it should immediately be re-surveyed with the same effort (Ojanen et al. 2013). Owing to this safeguard mechanism the majority of the patches were surveyed twice in 2018, reducing the potential influence of the ECE on nest detectability.
The documentation of the demographic crash of the Glanville fritillary metapopulation highlights the vulnerability of natural populations and underscores the importance of long-term monitoring programs to effectively document the consequences of ECEs on species, population and community dynamics. Our results further demonstrate that ECEs can impede conservation efforts by reducing the value of predictive models (see Matter & Roland 2017). For example, even within this extensively understood study system any effort to potentially mitigate the impacts of the ECE on local populations would likely have failed, since we would not have been able to predict which local populations were most vulnerable. The long-term implications of the ECE of 2018 are yet to be realized, but rapid recovery might be possible due to the high fecundity of the species (Saastamoinen 2007). Continuation of the metapopulation survey during the following years could therefore reveal whether and how the remaining populations bounce back from this rare disruption in extinction-colonization dynamics. This then also provides a unique and exciting     (see table S3). Meanwhile, these covariates captured extinction risk reasonably well in other years. Models were trained on ten random subsets of 80% of the available data. Plotted values correspond to mean values and line segments to standard error bars. As a general rule an AUC score of 0.7 is considered poor discrimination, 0.7 AUC score < 0.8 ≤ ≤ is considered fair, and 0.8 AUC score < 0.9 is considered good discrimination ≤ (Swets 1988