Mechanisms underpinning community stability along a latitudinal gradient: Insights from a niche‐based approach

At large scales, the mechanisms underpinning stability in natural communities may vary in importance due to changes in species composition, mean abundance, and species richness. Here we link species characteristics (niche positions) and community characteristics (richness and abundance) to evaluate the importance of stability mechanisms in 156 butterfly communities monitored across three European countries and spanning five bioclimatic regions. We construct niche‐based hierarchical structural Bayesian models to explain first differences in abundance, population stability, and species richness between the countries, and then explore how these factors impact community stability both directly and indirectly (via synchrony and population stability). Species richness was partially explained by the position of a site relative to the niches of the species pool, and species near the centre of their niche had higher average population stability. The differences in mean abundance, population stability, and species richness then influenced how much variation in community stability they explained across the countries. We found, using variance partitioning, that community stability in Finnish communities was most influenced by community abundance, whereas this aspect was unimportant in Spain with species synchrony explaining most variation; the UK was somewhat intermediate with both factors explaining variation. Across all countries, the diversity–stability relationship was indirect with species richness reducing synchrony which increased community stability, with no direct effects of species richness. Our results suggest that in natural communities, biogeographical variation observed in key drivers of stability, such as population abundance and species richness, leads to community stability being limited by different factors and that this can partially be explained due to the niche characteristics of the European butterfly assemblage.


| INTRODUC TI ON
A well-studied question in ecology is how species diversity (richness) influences the stability of the community (Elton, 1958;May, 1972;Tilman, 1999;Tilman & Downing, 1994). Community stability is known to be driven by two main factors: average population stability and synchrony (Doak et al., 1998;Thibaut & Connolly, 2013;Tilman et al., 1998;Yachi & Loreau, 1999). Briefly, average population stability describes the variability of each population over time and can be influenced by mean-variance scaling, where larger populations are relatively less variable (Kilpatrick & Ives, 2003;Taylor, 1961), as well as density-dependent (e.g. pathogen-induced mortality) and density-independent processes (e.g. mortality induced by extreme weather events). Population synchrony constitutes positively correlated population dynamics, which can be reduced by species interactions (e.g. competition) or through varying responses to environmental conditions. Average temporal synchrony between the populations in a given community is expected to decrease with species richness (Thibaut & Connolly, 2013), providing the main mechanism of the species-diversity relationship (Valencia et al., 2020), but associations between richness and average population stability are more variable (Jiang & Pu, 2009).
In addition to theoretical interest, the drivers affecting community stability are important to understand due to current concerns over biodiversity loss (IPBES, 2019) and because community stability ultimately underpins the stability of ecosystem function (Cardinale et al., 2012;Millennium Ecosystem Assessment, 2005). Butterflies are some of the best-studied insect groups and show evidence of declines in Europe (Warren et al., 2021) and North America (Forister et al., 2021). In the short term, reductions in community stability through lower species richness and/or mean abundance may result in fluctuations in ecosystem function (Greenwell et al., 2019) and over the longer-term reductions in mean functional provision (Weise et al., 2020). Consequently, identifying the drivers influencing community stability is likely to improve the effectiveness of environmental management.
Although theory around community stability is well developed, most testing has been conducted with plant or aquatic communities (Craven et al., 2018;Dıáz & Cabido, 2001; van der Plas, 2019).
For example, asynchrony increases with species richness in grassland communities (Isbell et al., 2019;Roscher et al., 2011), and asynchrony is the main driver of stability in plant communities across a range of habitats (Valencia et al., 2020). To date, there has been less work on natural animal communities, though asynchrony, along with additional effects of mean abundance and diversity, has been shown to be important to the stability of arthropod, bird, and bat communities across different habitat types in Germany (Blüthgen et al., 2016).
Similarly, population stability and asynchrony have a comparable effect size on the stability in bat, bird, and butterfly communities in France (Olivier et al., 2020), with both these findings consistent with results from plant communities.
A difficulty when considering communities at larger scales is that we might expect that some of the factors that influence population stability and asynchrony (e.g. abundance and species richness) may themselves vary substantially across space and time.
For example, the abundant-centre hypothesis posits that species will be most abundant at locations near the centre of their range (Andrewartha & Birch, 1954;Brown, 1984;Lawton, 1993), and so sites located near many species' range centres should have higher population stability due to mean-variance effects (i.e. larger populations being relatively less variable over time; Oliver et al., 2012Oliver et al., , 2014. Similarly, biogeographical factors such as latitudinal gradients (Hillebrand, 2004) or differences in isolation and land area (MacArthur & Wilson, 1967) can affect equilibrium species richness. average population stability. The differences in mean abundance, population stability, and species richness then influenced how much variation in community stability they explained across the countries. We found, using variance partitioning, that community stability in Finnish communities was most influenced by community abundance, whereas this aspect was unimportant in Spain with species synchrony explaining most variation; the UK was somewhat intermediate with both factors explaining variation.
Across all countries, the diversity-stability relationship was indirect with species richness reducing synchrony which increased community stability, with no direct effects of species richness. Our results suggest that in natural communities, biogeographical variation observed in key drivers of stability, such as population abundance and species richness, leads to community stability being limited by different factors and that this can partially be explained due to the niche characteristics of the European butterfly assemblage.

K E Y W O R D S
Bayesian analysis, biodiversity, butterflies, community stability, diversity-stability, niche, population stability, variance partitioning As well as changes in the mean levels of abundance and richness across different regions, we might also expect within-region variation depending on factors such as habitat structure and environmental heterogeneity. When evaluating observational data, this can lead to differences between the variance explained by a factor and its effect size, for example, even if species richness has a large effect on stability, it cannot explain observed variation in community stability in regions where species richness is constant across sites.
To tackle this complexity, we utilise the ecological niche concept (Hutchinson, 1957) to first understand how the underlying species characteristics of the European butterfly assemblage might influence population abundance and stability, and then look at how these factors, and their variability, subsequently affect community stability in butterfly communities across three countries (Spain, Finland, and the UK). Fundamental to niche modelling, at both the species and community scales (Hirzel & Le Lay, 2008;Poggiato et al., 2021), is that species occurrence under local conditions is dictated substantially by the n-dimensional environmental niche space. Consequently, notwithstanding nuances of extinction debt and colonisation credit (Kuussaari et al., 2009;, species richness should be predictable from the location of a site relative to the niches of the species in the regional species pool. Similarly, the mean abundance should be informed by site position relative to the niches of the species at the site (i.e. the abundance-niche-centre hypothesis; Martínez-Meyer et al., 2013;Yañez-Arenas et al., 2014) as has been found for birds (Osorio-Olvera et al., 2020). Species at the edge of their niche may also have lower average population stability (Mills et al., 2017;Oliver et al., 2012) and species with larger niche breadths may be more robust to local environmental variation. Finally, when a site is in a different niche position for two species, they may have differing responses to environmental variation at the site -reducing temporal synchrony in population dynamics (hereon referred to simply as synchrony). Therefore, the niche provides a way to link species characteristics and community dynamics. We apply this here to evaluate how the niches of the European butterfly assemblage impact community dynamics across three countries (Spain, the UK, and Finland) that span bioclimatic regions (Schmucki et al., 2016). Using structural equation modelling, we then evaluate the relationships between the proposed main mechanisms (community abundance, species richness, synchrony, population stability) and community stability.
Finally, we explore how much variation each of these factors explains across the countries for these naturally occurring populations.
A computationally fast approach is used to generate bioclimatic niche hypervolumes (Blonder et al., 2014(Blonder et al., , 2018 for each species and then niche-based metrics are derived to test potential factors affecting first species richness, mean abundance, and then the mechanisms of community stability. Hypothesised relationships ( Figure 2) are tested in a piecewise hierarchical Bayesian structural equation model that we designed to robustly account for between country variation in relationships and within country spatial autocorrelation. After fitting these models, we then use a Bayesian approach to variance partitioning (Schulz et al., 2021) to explore the importance of the different mechanisms for explaining community stability.

| Butterfly data
Data were collected with Butterfly Monitoring Schemes running in Finland, Northern Spain (Catalonia) and the UK. In these schemes, volunteers count butterflies along line transects following a standardised framework called the 'Pollard walk' (Pollard & Yates, 1993).
We processed the counts using generalized additive models (GAMs) to derive an index of abundance per site and year Schmucki et al., 2016). Note that a single number is provided both for multivoltine and univoltine species; therefore, it is the total abundance of adult butterflies across the observation period. If there are missing counts at a site, the GAM model interpolates abundance based on counts and phenological patterns observed at other sites in the same bioclimatic zone (Metzger et al., 2013), producing unbiased estimates of abundance (Schmucki et al., 2016).
In the different countries, the number of sites and the number of years of monitoring varies: Finland (1999:2017, number of sites = 107), Spain (1994Spain ( :2017, and the UK (1976:2017, n = 2128). Therefore, the country with the shortest running scheme, Finland, set the time period for our study (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017), and only sites with ≥10 years of data were retained, leaving 59 and 55 sites from Spain and Finland, respectively. To maintain a balance in the sampling between the countries, 60 sites with ≥10 years of data were randomly sampled from the larger collection of UK sites under the condition that each was >20 km from the others. To ensure robust estimates (restrict the risk of bias) from the GAM, annual indices of abundance with >50% of missing observations were also removed; this retained the same number of transect sites (n = 174) but removed some years at sites when the counts were made infrequently. A further 18 sites were removed prior to the final analysis, as they either contained only one species or too few occurrences per species to calculate average population stability (each species only seen once during observations); note this filtering occurred after some species were removed due to data limitations when constructing the hypervolume (see below). This left a final dataset of 156 sites (Spain = 54, Finland = 48, UK = 54).
We next utilised 19 'bioclimatic' variables provided by WorldClim 2.1 (Fick & Hijmans, 2017) at a resolution of 2.5 arc minutes and clipped the data to the approximate region of the species observations (i.e. western Europe; Figure 1). WorldClim's bioclimatic variables are derived from temperature and precipitation records and represent both mean and extreme conditions at a location. We fitted a principal components analysis (PCA) to the bioclimatic variables. The first six axes were retained as they represented >95% of the variance across the 19 bioclimatic variables. For each species' observations, we collected the first six coordinates projected onto this PCA space and used this to fit a 6-D hypervolume using multidimensional kernel density estimation (Blonder et al., 2014) using the 'hypervolume' 2.0.12 package (Blonder et al., 2014(Blonder et al., , 2018. In all, 12 species were removed from all further analyses at this point as there was an insufficient number of observations to fit the hypervolume algorithm (minimum 50 observations); thus, we retained data for 145 species in total. We expect the removed species, being very rare, to have little impact on overall community properties.

| Niche and community metrics
Five models (M1-M5) were used to link the niche metrics to community properties based on our a priori hypotheses. This required the calculation of four niche-based metrics and five measures of community properties associated with stability. For clarity, we provide a diagram of our hypotheses in Figure 2.
For M1, species richness is predicted to increase if a site is located nearer to the centre of the niches of more of the total species pool (i.e. all 145 species). Species richness was measured as the number of species observed at the monitoring site during the study period. To construct the average distance between the site location in niche space and the species pool niche centres, we took the average Euclidean distance between the site's location in a 6D hypervolume and the niche centroid of every species, hereafter referred to as overall niche distance. Thus, a site with a larger overall niche distance is further away from the niche centre of the total species pool.
M2 predicts that mean population abundance will increase as species are nearer the centre of their fundamental niche, and will increase for species with larger niches, that is environmental generalists are on average more abundant. To calculate the distance of species to their niche centres, we took the weighted average Euclidean distance between the niche centroids of each species found at a site and the location of the site in 6D niche space. The weighting means that the distance of the most abundant species in a community is weighted most highly; we term this measure niche mismatch. To calculate the size of species' environmental niches, we took the abundance-weighted mean volume of the hypervolumes which we term mean niche volume. Mean abundance was calculated from the density observed at each site. This metric was preferred to the original index as it accounts for differences in transect length between study sites. We log-transformed the mean abundance to better meet the assumptions of linear modelling. We also use both mean abundance and community mean abundance in our analysis as the community abundance might affect community stability (M5 see below), whereas mean abundance tests the link between niche-mismatch, weighted to account for evenness, and population abundance.
For M3, we predicted that average population stability will be highest when the site is closer to the centre of the species niches (niche mismatch), when the niche breadth of species is larger (mean niche volume), and when the abundance of the species is higher (mean abundance). To calculate average population stability, we follow Thibaut and Connolly (2013): Here m(i) refers to the mean abundance of species i and mc refers to the sum of species mean abundances in the community; thus, the score is abundance weighted. The second term is the mean of the species abundance divided by the standard deviation-the inverse coefficient of variation (Tilman, 1999).
In M4, we predict that synchrony will decrease with species richness and where species have lower niche-overlap. To measure syn- Map of the study extent, colours represent average annual temperatures (°C; bioclimatic variable 1) and black circles show butterfly monitoring sites, shaded to show any overlap of the points used to represent sites.
Here i and j refer to species i at site j, the numerator is the sum of all elements of the covariance matrix of the species at a site, and the denominator is the species-level variances in the presence of perfect synchrony. This measure is therefore standardised, accounting for differences in richness and variance. The synchrony index always takes a score between 0-no synchrony or population variance, and 1-perfect synchrony. Other measures of synchrony have been produced that have some advantages over φ in certain conditions (Gross et al., 2014); however, φ is easily interpreted and its result is often highly correlated with other methods.
To measure overlap between niches, we utilised the Jaccard Similarity which is the intersection of a pair of species niches divided by their union (Jaccard, 1912), once again always scoring between 0 and 1. This metric was calculated for all pairwise species comparisons and the mean was used to give a niche-overlap score.
In M5, we predict that synchrony and average population stability will be key to explaining differences in community stability, along with possible direct effects of species richness (Valencia et al., 2020) and mean community abundance (Blüthgen et al., 2016;Wagg et al., 2018). For mean community abundance, we took the mean of the summed density for all species across all years and we again log transformed to fit linear model assumptions. For community stability, we utilise the inverse coefficient of variation with a common correction for differences in sample size.
where μ refers to mean abundance for all species at the site, V to the variance, and n to the sample size. (2)

| Statistical analysis
To test the structure proposed in Figure 2, we produced five models (M1-M5) and connected them using a piecewise structural approach.
Each model was fit with a similar hierarchical model structure (see below).

Linear model:
Country-level parameters: Regional hyperparameters: Spatial offset: where i refers to site, l to other sites within the same country, and j to an index for each country. α refers to the intercept, β 1−n to coefficients of predictor variables 1 − n, X 1−n to data columns 1 − n, μ to mean, σ to standard deviation, and ρ is the correlation between slopes and intercept. For modelling spatial autocorrelation, η 2 refers to the maximal covariance between any sites within the country, p 2 to the rate of decline in correlation with distance, D il to the Euclidean distance between site i and l, and finally σ 2 i refers to within site (i = l) covariation with δ il set to 0 for i ≠ l and 1 for i = l.
Hierarchical or multi-level modelling is a probabilistic modelling framework where information between the different units (here the different countries) is pooled during estimation (Gelman et al., 2013;Gelman & Hill, 2006). Because the unit level parameters are drawn from the distributions of the population-level parameters (population in the statistical sense i.e. hyperparameters) which are estimated simultaneously from the data, this approach can lead to conservative estimates of effect sizes but reduces overfitting and therefore improves generalisability. The hierarchical model also takes account of the sampling structure, that is repeated sampling within units. Our model is a random slope and intercept model for each country as we assume that the mechanisms across the countries are largely similar, but there may be some differences depending on local conditions. In our models, all unit-level parameters are drawn from a multivariate normal across the hyperparameters which accounts for correlations between the slopes and intercepts. We also use a Gaussian process (Neal, 1998) to account for possible spatial autocorrelation in the dependent variable by fitting a spatial offset (K). The spatial offset for each site is drawn from a multivariate normal fitted to each country where the covariance matrix between sites is estimated using Euclidian distance and a squared exponential covariance function (L2 norm; McElreath, 2020a).
All variables were centred and scaled before fitting. We used weakly regularizing standard normal priors for the linear model parameters, an exponential prior with λ = 1 for the standard deviations, and η and p in the spatial offset. An LKF(2) prior (Lewandowski et al., 2009) was used on the correlation matrix for ρ.
To evaluate the dependence structure of the overall structural equation (Figure 2), we used directed separation tests (Shipley, 2000).
This consists of testing the implied independence between the variables given the structure presented in Figure 2. This is effectively a list of linear models (n = 23) where the independence between two variables is tested after controlling for other variables. Since our models were fitted within a Bayesian framework, we followed the approach of Clough (2012) and fitted each model and then approximated a p-value from the posterior of the parameter representing the dependence claim. p-values were derived from the posterior using the two-tailed test described in Shi and Yin (2021) and the independence of the whole structure was assessed from Fisher's C statistic (Fisher, 1954;Shipley, 2000). As can be applied in a standard piecewise structural approach (Lefcheck, 2016), we removed some of the independence claims (n = 5) and fitted them as correlated errors (i.e. measure their correlation) as these variables were derived from similar measures or may account for latent variables. The similar measures were mean population abundance-mean community abundance and niche mismatch-overall niche distance. To account for a possible latent effect of habitat quality species richness-mean abundance, species richness-mean community abundance were fitted as correlated errors. Finally, niche-overlap and species richness were also correlated as niche-overlap necessarily decreases when more species are added to a community; rather than fitting a separate model for this uninteresting result, we measured the correlation between the two metrics. Note that we fitted M4 with and without niche-overlap to assess whether this correlation influenced the fit between species richness and synchrony. Since it had little impact on the fit, we retained it to test the original hypothesis. We used a similar approach for the relationship between synchrony and community stability as the relationship was non-linear (results). In this case, we fit an alter-

| Variance partitioning
To measure the amount of variance in community stability explained by each of the independent variables in M5, we followed the approach presented in Schulz et al. (2021) for constructing variance partition estimates that take account of posterior uncertainty. For each country, we provide two measures of variance partitioning, first the variance partition, which is the proportion of variance in the dependent variable due to the linear term, and second the diagonal partition, which is the proportion of independent variation explained by the linear term. We use two measures because variance partition estimates can be biased by correlations in the independent variables and it is therefore recommended to provide more than one estimate (Schulz et al., 2021). We calculated these values for each country-

| RE SULTS
Throughout the results and discussion, we use evidence language  Figures S4 and S5).
For M1, there was evidence that the overall niche distance (the average distance between the site location in niche space and the species pool niche centres) of the community was associated with reduced species richness in Spain but little evidence for an effect in either the UK or Finland (Figure 3a; Figure S2). In M2, we evaluated factors influencing mean population abundance. We predicted that niche mismatch (the average distance of species to their niche centres) would decrease mean abundance and an increased mean niche volume (average niche hypervolume, indicating the degree of generalism, across species) might increase mean abundance. There was little evidence that niche-mismatch had any influence on mean abundance in any country and only in the UK was there evidence that mean niche volume increased abundance-but this effect also seemed to be driven by only a few sites ( Figure 3b) and should be interpreted with caution.
For M3, we predicted increased niche-mismatch would decrease population stability, and that mean abundance and mean niche volume would increase population stability. There was evidence for the effect of niche-mismatch in the UK and weak evidence for sites in Spain and Finland (Figure 3d). The estimated effect sizes (M3) were similar across all countries. Similarly, there was evidence that mean abundance increased population stability in the UK and weak evidence for Spain and Finland. Again, estimated effect sizes were similar across all countries (Figure 3c). There was little evidence that mean niche volume influenced population stability in any country ( Figure S3).
For M4, we predicted species richness will reduce synchrony and niche-overlap will increase it. There was evidence that species richness decreased synchrony in Spain and Finland and weak evidence for this effect in the UK (Figure 3e). Effect sizes were similar across all countries. There was little evidence that niche-overlap influenced synchrony in Spain and Finland. In the UK, contrary to expectation, there was weak evidence that niche-overlap decreased synchrony.
Nevertheless, this effect may be driven by only a few communities and should be interpreted with caution ( Figure S3).
Finally, for M5 we predicted that population stability, mean community abundance, and species richness would increase community stability and synchrony would decrease it. We found evidence that population stability increased stability in Spain and the UK, but little evidence that it affected the communities in Finland (Figure 3g).
Contrastingly, we found evidence that mean community abundance was associated with increased community stability in Finland and the UK, but little evidence that it affected community stability in Spain (Figure 3f; note, after accounting for the non-linear association between synchrony and community stability, the relationship between mean community abundance and community stability was estimated closer to zero for Spain). We found evidence that synchrony decreased community stability in all countries ( Figure 3h) and there was little evidence of a direct effect of species richness for any country ( Figure S3). At the regional level (i.e. statistical population, estimating a generalised result for all European butterfly communities), due to variation in responses between countries, there was much higher uncertainty and only the effects consistent across countries provided any evidence towards general responses ( Figure S1): for M3, there was weak evidence that niche-mismatch decreased average population stability and that mean abundance increased it. For M4, there was weak evidence that species richness decreased synchrony.
Finally, for M5, there was weak evidence that synchrony decreased community stability and mean community abundance increased it, whereas there was little evidence for an effect of average population stability.
Although not explicitly tested, we also noted differences in the range of some of the independent/dependent variables between countries. Finnish communities tended to have lower average population abundance (Figure 3c,f; intercept in Figure S2b) and mean community abundance and species had lower mean niche volumes and less niche-mismatch between the species and the sites. On the other hand, the UK had comparable abundance measures to Spain but the lowest species richness. Spain had communities with the highest species richness ( Figure S6) and species spanning the greatest range of niche space (Figure 3e; intercept in Figure S2a). Communities in the UK also had the highest niche-overlap and Spain the lowest ( Figure S3).

F I G U R E 3
Selected marginal effects from tests in M1-M5. Solid lines show means from the posterior and dotted lines 89% credible intervals (McElreath, 2020a), dashed lines are used instead of solid lines when there is little evidence of an effect, that is the 80% posterior credible interval contains zero; (a) Species richness against overall niche distance, (b) Log mean abundance against mean niche volume, (c) Average population stability against log mean abundance, (d) Average population stability against niche mismatch, (e) Synchrony against species richness, (f) Community stability against log mean community abundance, (g) Community stability against average population stability, (h) community stability against synchrony, (i) representation of the spatial covariance for each country (Spain, Finland, and the UK; using the same colours as for panels a-h), line shading represents the correlation between sites; note, countries are not to scale.
There were differences between the countries in the variance partitions for the independent variables in M5. In Spain, synchrony explained the largest amount of variance, followed by average population stability, while mean community abundance and species richness explained very little variation (Figure 5a). For Finland, mean community abundance explained most variation, though there was a large amount of uncertainty around its contribution, while synchrony explained a much smaller amount of variation, and average population stability and species richness explained very little variation (Figure 5b).
For the UK, aside from species richness, all variables explained some variation though synchrony tended to explain the most (Figure 5c).
The spatial offset accounted for a similar proportion of variation across countries, and the UK also had the largest amount of unexplained variation (this can also be observed in the fits in Figure 3f-h).

| DISCUSS ION
Using a structured hierarchical approach, we explored first how the niche characteristics of the European butterfly assemblage impacted key drivers of community stability (species richness, abundance, population stability, and synchrony) and how these, in turn,

F I G U R E 4
Overview of results in M1-M5 for (a) Spain, (b) Finland, and (c) the UK. Black lines represent positive and red lines negative associations with either evidence or weak evidence. Dashed lines in grey are used when there is little evidence for an effect, that is the 80% posterior credible interval contains zero. Correlated errors are shown in Figures S4 and S5. explained differences in community stability in butterfly communities in different regions of niche-space (namely Spain, Finland, and the UK). We found the niche metrics we constructed were somewhat informative as they influenced species richness and average population stability but not all our hypotheses were supported and results between countries often varied (Figures 3a-d and 4). The factors directly influencing community stability were more consistent, though there were differences between the countries with stronger evidence for an effect of mean community abundance on stability in Finland and the UK compared to Spain and little evidence for an effect of population stability in Finland. The amount of variance explained by these factors also differed between the countries.
Synchrony explained variation in community stability consistently across all three countries, though it explained a considerably larger proportion of its variance in Spain. In contrast, population stability explained little of the variation observed in Finnish communities but a consistent proportion in the UK and Spain, though less than synchrony in both cases. Finally, mean community abundance had a direct effect on community stability in Finland (here explaining most variation) and the UK but no effect in Spain. In all three countries, species richness showed no direct effect on community stability.
We discuss first the mechanisms operating at a biogeographical scale that impact species richness and then evaluate how species niches subsequently influenced specific stability mechanisms (M1-M3). There were large differences in species richness between the countries, Spain had higher average species richness and communities with an equivalent species richness to Finland and the UK were at locations less climatically suitable (higher overall niche mismatch) relative to the niches of European assemblage (Figure 3a).
The Mediterranean basin is recognised as a particular hotspot for insect diversity (Stefanescu et al., 2004) and the observed latitudinal gradient is also consistent with global patterns of butterfly species richness (Pinkert et al., 2022). This suggests that some of the main drivers of differences in community stability, via effects of species richness, are likely a combination of niche-based and historical factors. Two main reasons may explain the species richness gradient: post-glacial dispersal limitation and tropical niche-conservatism.
In the quaternary period, repeated glaciations caused many species to persist only in southern refugia (Hewitt, 1996); consequently, the lower species richness in Finland, the UK, and other northern European countries could be a result of historical climate and dispersal limitation (Essl et al., 2011), particularly for the UK as islands have increased risks of local extirpation and reduced immigration, limiting recolonisation success (Konvicka et al., 2006;MacArthur & Wilson, 1967;Whittaker & Fernández-Palacios, 2007). Alternatively, tropical niche-conservatism suggests that insect groups originating in the tropics require adaptation to factors like freezing that occur in cooler temperate climates (Wiens & Donoghue, 2004) limiting species distributions. The UK and Finland are at the northern range limits of many species owing to their cooler climates (Warren et al., 2001) and it is possible a single niche-dimension, such as temperature Melero et al., 2022) or evapotranspiration rates (Hawkins, 2010;Hawkins & Porter, 2003), may disproportionately influence persistence, particularly during extreme events. This would reduce species in the UK and Finland to an adaptable and resilient subset of the wider European assemblage (Thomas, 1993).
The differences in species richness had a larger impact on the proportion of variance explained by mechanisms of community stability rather than their effect size. Synchrony had the greatest effect and provided the strongest evidence for a generalised impact on community stability (though still weak evidence at the regional level). The effect was estimated to be largest in Spain and lowest in Finland; however, the non-linearity in this response and lower levels of community stability in Finland might mean the relationship between synchrony and community stability is more consistent than we observe, as expected from theory (Thibaut & Connolly, 2013;Wang & Loreau, 2014). Similarly, we found evidence for a positive effect of average population stability on community stability in the UK and Spain, but little evidence of an impact in Finland. The Finnish communities were clustered around lower values of average population stability (Figure 3g) possibly limiting our capacity to measure the effect, but the varying fits for the UK and Spain also suggest that the role of average population stability may be more variable than synchrony in butterfly communities.
Mean community abundance led to increases in community stability in Finland and the UK but there was little evidence of an impact in Spain. Once again, differences in the range of abundances may explain between country variability as the Finnish sites had the lowest abundances with only the most abundant sites overlapping with the lowest abundance of Spanish sites ( Figure 3f). As with population stability, in general, we expect higher stability from more abundant communities due to mean-variance scaling relationships, that is Taylor's power law (Kilpatrick & Ives, 2003;Taylor, 1961).
Alternatively, there may be some effect of a latent variable, such as habitat quality increasing both abundance and community stability. The weak evidence we find at the regional level (i.e. the entire statistical population) for mean community abundance and the little evidence for average population stability suggest the general effects of these factors remain uncertain. Finally, neither at the country nor the regional level do we find evidence that species richness influences community stability after accounting for its effect on synchrony. The evidence for the association we find between synchrony and species richness suggests that the diversity-stability effect in butterfly communities is mainly mediated through synchrony.
The amount of variation explained by the key mechanisms of community stability differed between countries. More variation was explained by synchrony in Spain than in the UK or Finland, which was likely driven by the greater variation in species richness in Spain ( Figure 3e). Part of this may be because the species pool in Spain is larger, but also that the Catalonian scheme in the north-eastern tip of Spain (Figure 3i) covers part of the Pyrenees with a large range of altitudes and varying amounts of woodland cover (García Viñas et al., 2006). In contrast, the largest amount of variation for Finland is explained by mean community abundance which, as we suggested, might be due to a latent effect of habitat quality or how average population stability varies within and between Finnish communities.
A consistent result we do find across all three countries is that synchrony explains more variation than average population stability and that species richness is explaining very little variation.
Our results regarding our niche-based hypotheses (M1-M3) were mixed. For Spain, we found evidence that species richness was higher if the site was located closer to the niche centre of the European assemblage. This result is expected as niche position is related to the probability of persistence and occurrence; a core principle of niche-distribution modelling (Elith & Leathwick, 2009). We note here that a simple aggregate score can be informative about the species richness across the community in some circumstances though the lack of evidence for it predicting richness in the UK and Finland suggests that more sophisticated approaches may be required to link community characteristics to the niches of the species pool.
We also anticipated, based on the abundant-centre hypothesis (Andrewartha & Birch, 1954;Brown, 1984;Lawton, 1993)  We did find evidence for an effect of niche-mismatch on average population stability. Populations at the margins of species niches being more variable is a natural extension of the abundant nichecentre hypotheses, and populations at geographical range edges have been shown previously to be more variable (Mills et al., 2017;Oliver et al., 2012). We also provide weak evidence that this effect can be generalised, suggesting consideration of species position within their climatic niche could be important for understanding population stability. This may be particularly so given the recent, and projected, increase in extreme climatic events (Donohue et al., 2016;Ummenhofer & Meehl, 2017). We also found evidence that mean abundance increases population stability, as expected from meanvariance relationships (Kilpatrick & Ives, 2003;Taylor, 1961), and again the consistency between countries meant there was some weak evidence supporting that this effect may be general.
We expected mean niche volume to increase average popula- would be predicted to have lower scores on a niche axis related to precipitation or aridity. Finally, niches are constructed with uncertainty and so including more data on occurrence may lead to better estimates of niche metrics, particularly for rarer species.
In summary, across three European countries, we have explored niche-based factors that influence the key drivers of community stability (community abundance, species richness, synchrony, population stability), and then evaluated how these, in turn, influenced community stability. We found that ultimately, community stability is affected by broad biogeographical factors influencing the richness of the species pool across our countries, but that niche characteristics can modify the resultant stability of the community. There was some variability in the effect size of the drivers of stability across countries, but large differences in the proportion of variation they explained. Synchrony was the most consistent mechanism explaining differences in community stability and accounted for the diversity-stability relationship in European butterfly communities.
Regarding environmental management, our results suggest that efficiently targeting improvements to community and ecosystem function stability is likely context specific and that managers could consider the achievable increases in a driver of community stability (e.g. average population stability) alongside its effect size. Finally, our results suggest considerations of niches can be informative for understanding community dynamics. At large scales, the underlying differences in mean and variability of the factors that constitute the building blocks of community stability (e.g. species richness and abundance) mean their importance for community stability varies and that stability can occur through alternative processes dependent on regional context.

AUTH O R CO NTR I B UTI O N S
Luke Christopher Evans designed the study and performed the analysis, Reto Schmucki processed the data prior to analysis, Luke Christopher Evans, Tom Henry Oliver, and Yolanda Melero wrote the first draft, and all other authors contributed substantially to revisions.

ACK N OWLED G M ENTS
We thank the volunteers collecting butterfly data and the funders of the schemes for obtaining the data required for this study.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflicts of interest.