Bioclimatic context of species' populations determines community stability

Evans, L. C. ORCID: https://orcid.org/0000-0001-8649-0589, Melero, Y. ORCID: https://orcid.org/0000-0002-4337-1448, Schmucki, R., Boersch Supan, P. H. ORCID: ‐ https://orcid.org/0000-0001-6723-6833, Brotons, L., Fontaine, C., Jiguet, F., Kuussaari, M., Massimino, D., Robinson, R. A., Roy, D. B., Schweiger, O., Settele, J., Stefanescu, C., van Turnhout, C. A. M. and Oliver, T. H. (2022) Bioclimatic context of species' populations determines community stability. Global Ecology and Biogeography. ISSN 1466-8238 doi: https://doi.org/10.1111/geb.13527 Available at https://centaur.reading.ac.uk/105396/


| INTRODUC TI ON
The relationships between diversity, stability and ecosystem functioning have long been debated (Elton, 1958;May, 1972;Tilman & Downing, 1994), and they remain in focus owing to recently observed (Díaz et al., 2019) and projected (Ceballos et al., 2020) declines in biodiversity. Understanding the responses of species and communities to environmental disturbance is also a priority for conservation, given the projected increases in the magnitude and frequency of extreme climatic events (Donohue et al., 2016;Ummenhofer & Meehl, 2017).
Currently, there are large concerns around insect declines (Cardoso et al., 2020), although debate about their severity and extent is ongoing (Simmons et al., 2019;Thomas et al., 2019). Butterflies, one of the most monitored and studied insect groups, show evidence of declines across Europe (Warren et al., 2021) and in North America (Forister et al., 2021). Consequently, understanding the drivers influencing the dynamics of populations and communities of this insect group is important, because synchronous fluctuations in abundance can lead to short-term losses of ecosystem function (Greenwell et al., 2019) and long-term declines to reduction in mean function provision.
Community and population stability can vary spatially. Butterflies, along with declines in abundance, exhibit recent distributional changes (Warren et al., 2021). Consistent with the Grinnellian niche concept (Grinnell, 1917a, Grinnell, 1917b, distributions of butterfly species are driven largely by abiotic factors, such as climate (Settele et al., 2008), and their population dynamics are driven by weather (Palmer et al., 2017;WallisDeVries et al., 2011). In combination, we can expect that populations are most abundant and stable near the centre of their niche range and are most sensitive to environmental variation at range edges (Brown, 1984;Mills et al., 2017;Osorio-Olvera et al., 2020).
Consequently, the local dynamics of communities might be impacted by the larger-scale biogeographical context of the constituent species, in addition to local factors, such as weather anomalies.
The overall stability of a community is contingent on several factors. Community stability often increases with biodiversity (Cardinale et al., 2012;Jiang & Pu, 2009), although the strength of this relationship can vary between systems (Campbell et al., 2011).
Early mechanistic explanations for diversity-stability relationships included the insurance hypothesis (Yachi & Loreau, 1999) and the portfolio effect (Doak et al., 1998;Tilman et al., 1998). These are now combined in modelling frameworks demonstrating that higher species richness increases asynchrony in population dynamics and, along with mean population stability, influences the overall community stability (Thibaut & Connolly, 2013;Wang & Loreau, 2014).
Therefore, when considering butterflies and other insect groups, it is necessary to consider factors that generate asynchrony or that lead to higher overall population stability. Asynchrony is often associated with competition, but this might be less important in mobile animal communities where competition is diffuse (Loreau & De Mazancourt, 2008). Many animal communities at similar trophic levels are characterized by low competition owing to limiting similarity (Macarthur & Levins, 1967) or specialization (Futuyma & Moreno, 1988). Butterflies are resource generalists as adults (Dennis, 1992;Sharp et al., 1974) but more often specialists as larvae, owing to the adaptions required to circumvent structural or chemical plant defences (Ehrlich & Raven, 1964), which might limit interspecific competition. Although interspecific competition is present in phytophagous insects (Kaplan & Denno, 2007) and can be present in some butterfly species with shared host plants (Millan et al., 2013), or through apparent competition (Audusseau et al., 2021), asynchrony might be driven primarily by differential responses to weather anomalies, given the known effects of climate and weather on distributions and population dynamics.
Likewise, weather can have large impacts on population stability, particularly in association with range position. The links between distribution, abundance and geographical range have long been of interest in macroecological theory (e.g., abundant centre hypothesis; Brown, 1984), often supported by mixed evidence (Sagarin & Gaines, 2002). But substitution of the geographical position with niche position (abundant niche centre hypothesis; Martínez-Meyer et al., 2013;Yañez-Arenas et al., 2014) shows support for increased abundance (Osorio-Olvera et al., 2020), positive population trends (Manthey et al., 2015) and genetic variation (de Mazancourt et al., 2014)  Main conclusions: Future extreme heat events will be likely to impact species negatively across their thermal range, but might be particularly impactful on populations at the hottest end of the thermal range. Thermal range position influenced community stability because range edge communities were stable. However, the prediction of community stability from thermal range position is challenging because of nonlinear responses to temperature, with small temperature anomalies producing weak compensatory dynamics, but large extreme events synchronizing dynamics.

K E Y W O R D S
asynchrony, biodiversity, biogeography, community stability, diversity-stability, insects, integrated Laplace approximation, long-term monitoring, range position range edges (Mills et al., 2017), although the effect of niche position is less explored. Consequently, the mean population stability of a community might be impacted not only by the absolute size of a weather anomaly, but also by the niche position of the population. Therefore, for both the main factors influencing community stability (i.e., asynchrony and population stability), the bioclimatic context of species and their responses to key weather variables might be key factors to understand and predict differences in community stability at large scales.
A challenge with selecting a key environmental driver for community stability is that species vary in their sensitivity to different environmental variables (Lawson et al., 2015). Functional responses to many variables, although identified (McDermott Long et al., 2017;Roy et al., 2001), are not typically known in detail, limiting our understanding of how responses combine to influence the whole community. Therefore, the approach applied here is to select a single influential type of anomaly and evaluate how species responses contribute to community stability, asynchrony, mean population stability or any other stabilizing mechanism. Likewise, to understand the impact of bioclimatic context, we apply the abiotic niche (Grinnell, 1917b;Hutchinson, 1957) and the abundant niche centre (Martínez-Meyer et al., 2013;Yañez-Arenas et al., 2014) concepts along a single major axis. This allows us to understand how the range position of species combines with local climatic anomalies to influence community stability at large scales. We select here temperature, owing to its influence on key biological rates that impact fitness (Kingsolver, 2009) and because it gives insight into the potential ecological impacts of climate change (Altwegg et al., 2017;Palmer et al., 2017;Pandori & Sorte, 2019).
Temperature responses are relatively well characterized in ectotherms (Angilletta et al., 2002(Angilletta et al., , 2010, with biological rates typically increasing up to an optimum temperature before subsequently declining rapidly (Briere et al., 1999;Huey & Kingsolver, 1989;Shi & Ge, 2010). Therefore, an expectation is that populations of species at the cold edge of their temperature range (position in the temperature range is termed here "thermal range position") will have contrasting responses to local temperature anomalies compared with populations of species that are at the hot edge (i.e., position in the thermal range is a potential mechanism driving asynchrony in population dynamics of different species; Jiguet et al., 2010). Position in the thermal range can also impact the average stability of a population because those near their optimum might have milder responses to anomalies than those at range edges where populations are prone to crashes (Mills et al., 2017;Oliver et al., 2014). However, temperature F I G U R E 1 Hypothesized responses to temperature anomalies at five sites across a thermal range for species showing different combinations of thermal specialism and local adaptation. The mean temperature at the five sites is indicated by the vertical dashed lines, with colour representing the mean temperature (e.g., blue to red, cold to hot mean temperatures). In ectotherms, thermal performance (e.g., rates of growth, development and survival) are characterized by a shape that increases to a maximum before decreasing more sharply. We define local adaptation as the matching of the thermal performance curve for a species to the mean local conditions. (a,c) Species with low local adaptation have a single curve across their range. (e,g) In contrast, adapted species have their maximum thermal performance at the mean conditions at a site. To generate the potential impacts of temperature anomalies on population growth rate for these adaptationspecialism combinations, we simulated temperature observations for the five sites: one site at the optimum temperature for the performance curve shown in (a,c), then at 80, 90, 110 and 120% of this value. Linear models predicting logarithmic growth rate were then fitted to these data with fixed effects matching the structure of a model described below in Equation (1). (b,d,f,h) The expected mean responses. Further details about the simulation are presented in the Supporting Information (Appendix S2) responses might also be influenced by local adaptation or acclimatization (Huey & Kingsolver, 1989), impacting the expected mean population stability and asynchrony in the community.
Here, we test how the position in the thermal range and temperature anomalies combine to influence community stability at a European scale. We test the following four hypotheses (H 1 -H 4 ): H 1 , populations at the cold and hot range edges of their thermal range will have contrasting responses to local temperature anomalies (see Figure 1 for a simulation articulating H 1 ); and H 2 , communities with a more mixed composition of species' thermal range positions will have more stable dynamics overall. The two further hypotheses test possible mechanisms driving the effects of thermal range position: H 3 , communities composed of species from a mix of range positions will have contrasting population responses to temperature anomalies, leading to asynchrony in abundance and higher overall community stability; and H 4 , populations near the centre of species' thermal ranges will be more stable, meaning that communities with more populations at the thermal range centres will, overall, be more stable (for a diagram articulating H 2 , H 3 and H 4 , see Supporting Information Figure S1). We used the inverse of temporal variation in total community abundance as a proxy for community stability (Tilman, 1999) because it is tightly related to ecosystem function owing to mass ratio effects (Dangles & Malmqvist, 2004;Grime, 1998;Smith & Knapp, 2003).

| Data
Data were collected with butterfly monitoring schemes carried out in Finland, Spain (Catalonia) and the UK. The schemes consist of a network of sites where volunteers count butterflies along transects following a standardized framework called the "Pollard walk" (Pollard & Yates, 1993). We then processed counts using generalized additive models (GAMs) to provide an index of abundance per site and year (Dennis et al., 2013;Schmucki et al., 2016). In the case of missing counts, we interpolated from the GAM fitted to counts made at other sites in the same bioclimatic zone (Metzger et al., 2013;Schmucki et al., 2016). This approach generates unbiased estimates of abundance and performs better than interpolations from simple linear regressions (Dennis et al., 2013). However, to assure robust estimates, we removed indices of abundance with >50% of missing observations (Schmucki et al., 2016). of sites = 2128). Therefore, Finland, the country with the shortest scheme, set the study period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). From each of the schemes, only those sites with >10 years of data were retained, leaving 59, 55 and 762 sites from Spain, Finland and the UK, respectively. The data included 157 species, and Spain had the highest species richness and the highest number of species (S) unique to the country (S = 131, unique = 74), followed by Finland (S = 58, unique = 12) and the UK (S = 58, unique = 8). We converted counts into densities by dividing the indices of abundance by the transect length, thus standardizing the measure of abundance across transects of different lengths.
Temperature data were obtained from the European Climatic and Assessment Dataset (ECAD) project (Haylock et al., 2008;Klok & Klein Tank, 2009). ECAD provides gridded daily temperatures at a 0.1° scale through interpolations from observations collected by a network of meteorological stations. The temperature was collected to the nearest degree at each site from the period of 1999-2017.

| Thermal range construction
The locations of the transect sites within the thermal range of each species were measured in environmental (temperature) space. We took this approach because our focus was the position of a species in its temperature niche, rather than geographical position, because changes in landscape features and seasonal variation in climate at a continental scale can decouple the relationship between geographical space and weather (Loarie et al., 2009).
It was also necessary to select the temperature for a given time period because, given our broad biogeographical scale, different locations will have different positions in the temperature niche over time (e.g., sites in Finland are warmer in summer than sites in the UK but colder in winter). Furthermore, given that our analysis included many species with varying phenologies, it required a broad period.
We selected spring (March-May) because it coincides with the developmental period for many species and should be a broad enough period to capture differences across our spatial scale without averaging out anomalies. To standardize the position in niche space, we first took the daily temperatures at each site for the spring period and calculated the mean for the preceding 30 years  to obtain an average site temperature (T). We then transformed these mean values (Equation 1) such that the maximum mean temperature for each species across all sites was given a score of one and the minimum mean temperature minus one. Thus, every site-by-species combination had a position within the thermal range between minus one and one: where T ij is the mean temperature of the site for species i at site j, and T is the set of mean temperatures from all sites for that species.
The thermal range preserved the relative difference between the temperatures at the sites where a species was observed but standardized thermal range differences between species, hence a score of +.75 covers 87.5% of the total projected temperature range for a species independent of the absolute size of the temperature range. Note that the temperature range across all our monitoring sites from Catalonia (Spain) to Finland will not always capture the full thermal range for every species (e.g., some extend beyond the southern Mediterranean region into Africa), but the standardization reflects the relative thermal range position in our sampled dataset.

| Hypothesis 1: Negative covariance in response to temperature anomalies at thermal range edges
We predicted that the population responses of a species, in terms of change in the logarithm of abundance between one year (N t−1 ) and the next (N t ), ln(N t /N t−1 ), to a local temperature anomaly will be dependent on the position of the site within the thermal range for the species and that populations at the cold and hot edges will have contrasting responses to a given direction of temperature anomaly (see articulation of the hypothesis in Figure 1). To test this hypothesis, we fitted the relationship between the change in population size in response to a local temperature anomaly (the difference each year between the observed seasonal mean temperature and the average seasonal temperature conditions at a site over the last 30 years) and the thermal range position for all species in a hierarchical Bayesian model.
We predicted, on average, a nonlinear response, with large temperature anomalies either above or below the mean reducing population size, but also that the shape/direction of this response might be influenced by thermal range position and would vary between species. Given that thermal performance might be nonlinear To start, we fitted the a priori plausible maximal model (Equation 2). Fitting this model produced singular fits in higherorder random effect terms. Consequently, we simplified the model by progressively removing the singular highest-order terms one at a time until we obtained a model with sensible precision overall random effects (Equation 2) ( Barr et al., 2013). The removal of these higher-order terms reduced the deviance information criterion (DIC; Spiegelhalter et al., 2002) and the Watanabe-Akaike information criterion (WAIC; Watanabe & Opper, 2010) scores of the model fit.
Here, subscripts refer to site j, species i and year k; σ 2 refers to the variance and f(θ xy ) the value of the Gaussian spatial field at location x,y; trange refers to the thermal range position and temp anomaly the local temperature anomaly.
We used the integrated Laplace approximation (INLA) method for approximate Bayesian inference (Rue et al., 2009), which provides a method for fitting a Matérn covariance spatial autocorrelation function, through a weak solution to a stochastic partial differential equation (SPDE) (Gómez-Rubio, 2020; Lindgren et al., 2011;Lindgren & Rue, 2015). Further information about prior selection is provided in the Supporting Information (Appendix S3).

| Hypothesis 2: Community stability from thermal range position
We predicted that communities with central thermal range positions on average will have more stable community dynamics. To measure community stability, we used the inverse coefficient of variation (CV) for yearly abundance counts summed across all species at a site across the sampling period (maximum of 30 years), hence it was a measure of variation in total community abundance. This approach is similar to stability measures used for studies of ecosystem functioning (e.g., Donohue et al., 2016;Hautier et al., 2015), although here we do not extrapolate our results to function. Given that time series were not identical in length, we used a bias correction for smaller sample sizes, because the CV can be underestimated from shorter time series, and we also limited the data to sites with ≥10 years of data (Equations 3.1 and 3.2): where † refers to the bias-corrected estimate of the coefficient of variation at site j and n to the sample size, σ to the standard deviation of abundance, and μ to the mean abundance.
Community stability is also known to be affected by species richness (Tilman & Downing, 1994) and mean abundance (Taylor, 1961).
Therefore, at each site, we calculated the mean thermal range position of the species (using range position from average yearly temperatures), the mean community abundance (measured as density to compare across transects of different lengths) and the species richness.
For the model, we predicted a quadratic effect of thermal range position on community stability, because stability should increase y ijk = lnN t−1 + temp anomaly + temp anomaly 2 + temp anomaly 3 + trange pos i + trange i : temp anomaly + trange i : temp anomaly 2 + trange i : temp anomaly 3 + lnN t−1i + temp anomaly i + temp anomaly 2 i + trange i : temp anomaly i + trange i : temp anomaly 2 with a mean range score of zero and decrease when species are, on average, at either side of the range edges, and we also included linear effects of species richness and mean community abundance to control for their effects on stability: Symbols used here are the same as in Equation 2.

| Hypothesis 3: Asynchrony
We hypothesized two mechanisms that might explain the effect of thermal range position on community stability: H 3 , communities composed of species from a mix of range positions will have contrasting population responses to temperature anomalies, leading to asynchrony in abundance and higher overall community stability; and H 4 , populations near the centre of species' thermal ranges will be more stable, meaning that communities with more populations at the thermal range centres will, overall, be more stable.
We predicted that if asynchrony was influenced by opposing thermal range positions, the mean synchrony in the community (i.e., averaged across all pairs of species) should be lowest in communities with central average range positions (i.e., concave up relationship). To test this, we calculated the average population synchrony (measured as the mean of all pairwise correlations across species' time series at a site), including only species pairs with ≥10 concurrent observations. We then constructed a linear model (Equation 5) for synchrony against mean range position and mean range position 2 . Given that synchrony is expected to increase with species richness (Ives et al., 1999;Thibaut & Connolly, 2013), this was calculated for each site and included in the model.
Here, μ indicates the mean, and other symbols are as in Equation 2.
Given that we were modelling an estimate of mean synchrony, the regression was weighted by the inverse of the variance to control for variation in the precision of the estimate of the mean owing to differences in sample size. INLA does not offer a straightforward method for implementing weighted regression;consequently, we used the "brms" package (Bürkner, 2017) to fit a weighted regression using Hamiltonian Monte Carlo through Stan (Carpenter et al., 2017).

| Hypothesis 4: Population stability
We predicted that average population stability might increase at the centre of the thermal range, leading to increased community stability. We calculated the stability of the population dynamics of each species at each site, as in Equations 3.1 and 3.2, along with the mean population abundance. Observations included only species-site combinations with ≥10 observations in total, as in Equation 4. To test our hypothesis, we constructed a model with fixed and random effects of range position, range position 2 and mean population abundance, along with random intercepts of year and site. After fitting, the year intercept was found to be near singular and was removed from the model. Population abundance was also Ln-transformed before the fit because there were a few large populations that had an outsized influence on coefficient estimates. Information on priors is provided in the Supporting Information (Appendix S1).
The hierarchical Bayesian models were fitted using the "INLA" package (Lindgren & Rue, 2015), and the package "brinla" was used to aid in diagnostics (Faraway et al., 2021). Polygons for barrier models were taken from the "rnaturalearth" package (South, 2017), and spatial fields were plotted using functions from the paper by Krainski et al. (2018). Code supporting the results is archived at DOI: 10.5281/zenodo.6350070.

| Hypothesis 1: Responses to temperature anomalies across the thermal range
The model fit did not conform precisely to a thermal performance curve, although responses of population growth to anomalies were nonlinear. Species in the warmer half of their thermal range were most impacted by large anomalies, showing decreased growth rates with hot anomalies and increased growth with cold anomalies (Figure 2a). The response to temperature anomalies was weaker overall at the cold edge of the range, and there was evidence of contrasting responses, with species at the hot edge performing worse with high-temperature anomalies and better with low-temperature anomalies compared with species at the cold edge (Table 1; Figure 2a). However, at extreme hot temperatures the populations from all range positions declined, and towards the warm edge the declines were substantial. The variance of the random effects (Table 1) suggested that there was interspecies variety in both the shape of the response and the importance of the interaction terms.
(6) Population stability ijk ∼ N y ijk | σ 2 , The spatial field demonstrated relatively broad-scale regional correlations (range c. 95 km) in average population changes in the areas such as the south-east of England, showing lower growth rates relative to much of the interior, and the north-eastern tip of Spain, showing higher growth rates than the interior (Figure 2b).  Table S1).

| Hypothesis 3: Asynchrony
We found that synchrony tended towards being slightly concave up with thermal range position (Figure 4a). However, species richness had a clearer effect, with species richness decreasing synchrony (Figure 4b,c; Supporting Information Table S2).

| Hypothesis 4: Population stability
For mean population stability, rather than an increase in the centre of the range, stability was lowest at the hottest range edge ( Figure 5a; Supporting Information Table S2), although the overall effect was weak. However, mean population abundance increased stability (Figure 5b).

| DISCUSS ION
We tested four hypotheses about how the position of species within their thermal range and anomalies interact to influence population responses and community stability: H 1 , species will have contrasting responses to temperature anomalies at different ends of their thermal ranges; H 2 , communities consisting of populations with central thermal range positions on average will be more stable; H 3 , communities composed of species from a mix of range positions will have contrasting population responses to temperature anomalies, leading to asynchrony in abundance and higher overall community stability; and H 4 , populations near the centre of species' thermal ranges will F I G U R E 2 (a) Marginal fit of the relationship of interannual population change [ln(N t /N t−1 )] to local temperature anomaly for different thermal range positions across all species. For the position in the standardized range shown on the colour bar, minus one refers to the coldest site occupied by the species and one to hottest site (95% confidence intervals for the fit are shown in Supporting Information Figure  S7). (b) Spatial field showing the change in intercept for interannual population change across the study area after accounting for the other effects (i.e., the field indicates areas with correlated dynamics and with higher or lower average growth rates) be more stable, meaning that communities with more populations at the thermal range centres will, overall, be more stable. Hypothesis 1 was supported in part, because we found that the location of species within the thermal range influenced responses to local temperature anomalies ( Table 1; Figure 2), and populations at the hot edge performed worse with high-temperature anomalies and vice versa.
However, responses were highly nonlinear relative to expectations (Figures 1 and 2), and with extreme anomalies the populations at the hot edge showed large declines (~60% decline for +3°C) with hot-temperature anomalies and population growth with cold anomalies (~30% increase for −2.5°C). Populations at the hot edge were also most responsive to temperature anomalies, with the cold edge showing smaller reductions in performance at both edges of the thermal range. Hypothesis 2 was also supported because, after accounting for species richness and community abundance, community stability decreased towards range edges (Figure 3; Supporting Information Table S1). The mechanisms (H 3 and H 4 ), however, were less clear, because synchrony decreased only marginally towards the centre, and mean population stability declined towards the hotter range edge but did not noticeably peak at the range centre.
Our analysis of responses to temperature anomalies (Figure 2a) suggests that the composition of thermal range positions at a site is unlikely to aggregate in a simple, consistent way to impact community stability. This is because the degree of synchrony between species responses might be contingent upon the relative size of the temperature anomaly. For example, two species at either side of the centre of their thermal ranges might have small asynchronous responses to moderate weather variation, but in an extreme hightemperature event, both species might be pushed beyond thermal limits (Sunday et al., 2014) and crash synchronously. Our test of the asynchrony mechanism supports this, because synchrony decreased only slightly in the centre of the range, with species richness having a much larger effect.
A few reasons might explain why we found only a weak impact of asynchrony with thermal range position. First, although temperature is an important driver, butterfly species are impacted by other weather variables, such as precipitation (Herrando et al., 2019;Roy et al., 2001) and aridity (Oliver et al., 2015). This might disrupt the expected opposing responses of species at different thermal range edges. Second, temperature can also affect species in unpredictable ways owing to its impact on natural enemies or host plants, and phenology might also vary across the study region, affecting the life stage impacted by the anomaly (McDermott Long et al., 2017).
Consequently, asynchrony might be a larger driver of stability in butterfly communities in general, but the contrasting responses are not strongly connected to the thermal range constructed here. Instead, large temperature anomalies might lead to synchrony rather than asynchrony.
Our results suggest that populations at the hot edge of the thermal range are, in general, more responsive to weather anomalies than those at the cold edge. Previous work has shown that butterfly populations are more variable at range edges (Mills et al., 2017;Oliver et al., 2012), but greater declines in population growth rates towards the hottest edge have also been noted in US butterflies (Breed et al., 2013) and birds (Jiguet et al., 2010). We provide further support for population instability at the hot edge with our ulation trends have been noted in many taxa, including butterflies (Martay et al., 2017); therefore, management might be required to enhance the stability of populations at the warm edge with future climate change (Oliver et al., 2015). Furthermore, the results suggest that relative niche position might be a simple but important indicator of the populations and communities most at risk from climate change (e.g., Settele et al., 2008).
Although populations towards the hot end of the thermal range were more unstable, the communities at these locations were only marginally less stable than in the centre, and communities in Finland, at the cold edge, were overall the least stable. This is likely to be attributable to a combination of these communities having lower abundance and species richness than the communities in the UK and Spain, respectively; Spain, in particular, has much higher species richness than the other two countries. We do not include asynchrony directly in Equation 4, which ultimately drives the species richness effect (Thibaut & Connolly, 2013), and we only include species richness, which could lead to an underestimation of the negative impacts of communities dominated by populations at the hot range edge. Therefore, a more structural account of the initial impact of species richness on asynchrony, then of asynchrony on community stability (e.g., Olivier et al., 2020), could be informative for understanding community stability at this broad spatial scale.
Furthermore, the compositional change in species between the countries might also account for some differences in community stability. The UK contains the northern range limits of many species owing to its cool climate (Warren et al., 2001), and the species assemblage that survives in the UK might be an adaptable and resilient subset of the wider European assemblage (Thomas, 1993).
Average trait differences between the assemblages might also ex- , effect sizes are shown as black squares with 95% confidence intervals as lines. An asterisk is used to highlight effects where confidence intervals do not overlap zero sensitivity to climate change for many species (Pacifici et al., 2017).
Finally, Spain also contains the largest number of unique species, and although there are large temperature differentials between sites in Spain, some species might have thermal range limits beyond those encountered. Consequently, extending the study to include a larger area and to include data from other European butterfly monitoring schemes might provide an increased resolution for understanding population and community responses to temperature anomalies at the European scale.
After accounting for the effects of thermal range position and temperature anomalies, we found differing patterns in the spatial fields that represent average rates of population change and community stability across the study area (Figures 2b and 3b). There is a general expectation that sites closer to one another should have similar dynamics because they share various environmental conditions (Tobler, 1970), which we found. The spatial pattern for rates of change of interannual population was in the range of 100 km, suggesting similar patches of average population change at regional scales for our countries. The spatial field for community stability, however, showed a smaller-scale patterning (range c. 25 km). The different spatial patterns are likely to be influenced by different drivers.
Population dynamics might be driven by regional weather conditions (Breed et al., 2013), possibly interacting with broad land-cover types (Stefanescu et al., 2011)  influenced by the cumulative effects of temperature (e.g., growing degree days), and in these cases the average temperatures or anomalies can be misleading (Denny, 2017;Sunday et al., 2014).
Linking community patterns with detailed mechanisms, such as context-dependent growth rates, is a major challenge and currently has been attempted only at small scales (e.g., White et al., 2020). In future, modelling efforts might be better able to connect species physiology to large-scale patterns observed for butterflies by using standardized schemes , but this would require tighter integration of the cycle of data collection, theory, modelling and testing than is typically practised (Boult & Evans, 2021;Dietze et al., 2018). Finally, our data might contain observational errors, which have been shown to explain variation in community stability (de Mazancourt et al., 2014). There has been exploratory work on the possible impacts of species detectability in butterfly counts (Isaac et al., 2011), but it remains unclear whether adjusting for detectability would improve estimates of real abundance. Given that we measure variation in relative abundance across years at fixed locations, where detectability should be reasonably constant, it is unlikely to have a large impact on our results.

F I G U R E 5
Marginal relationships between population stability and (a) thermal range position, or (b) natural logarithm of mean abundance. Spanish sites are shown in red, Finnish sites in blue and UK sites in green. Lines show marginal fits with 95% credible intervals. Note that (a) does not show the full extent of the residuals, in order to show the model fit better

| CON CLUS IONS
In summary, we have demonstrated that the interaction of thermal range position and local anomalies influences both interannual population change and community stability in butterflies at a broad biogeographical scale. We found that range edge communities were less stable and that populations near the hot edge of the thermal range were most strongly influenced by temperature anomalies. Responses were also nonlinear, meaning that many species might be impacted strongly by extreme heat events. Aggregating population responses across a single-niche dimension, although informative, does not simply predict community stability, because we found that small temperature anomalies can produce weak compensatory dynamics, but large extreme events might synchronize dynamics. Our results suggest that niche position is an important determinant of community and population stability, but the larger sensitivity of populations at the hot edge suggests that community stability at hot locations might be most impacted by climate change.

ACK N OWLED G M ENTS
We thank the volunteers for collecting their butterfly observations and the funders of the schemes for obtaining the data required for this study. The UK butterfly monitoring scheme is organized and funded by Butterfly Conservation, the Centre for Ecology and Hydrology, British Trust for Ornithology and the Joint Nature Conservation Committee. The Catalan butterfly monitoring scheme is funded by the Catalan Government, the Barcelona Provincial Council and other local partners. The Catalan butterfly monitoring scheme also incorporates the Andora butterfly monitoring scheme, which is run by Centre d'Estudis de la Neu i la Mutanya d'Andorra (CENMA) and funded by Govern d'Andorra. The Finnish butterfly monitoring scheme is organized and funded by the Finnish Environment Institute (SYKE). We also thank two anonymous reviewers and the handling editor Dr Petr Keil, who provided helpful comments for improving the manuscript.