Forecasts of butterfly future richness change in the southwest Mediterranean. The role of sampling effort and non-climatic variables

We estimated the potential impact of Global Warming on the species richness of Iberian butterflies. First, we determined the grid size that maximized the balance between geographic resolution, area coverage and environmental representativeness. Contemporary richness was modelled in several alternative ways that differed in how sampling effort was controlled for, and in whether the non-climatic variables (physiography, lithology, position) were incorporated. The results were extrapolated to four WorldClim scenarios. Richness loss is to be expected for at least 70% of the area, with forecasts from the combined models being only slightly more optimistic than those from the purely climatic ones. Overall, the most intense losses are predicted for areas of highest contemporary species richness, while the potential slightly positive or nearly neutral changes would most often concentrate in cells of low to moderate present richness. The environmental determinants of richness might not be uniform across the geographical range of sampling effort, suggesting the need of additional data from the least intensively surveyed areas. Re-assessing richness and its environmental determinants in the area proves necessary for more detailed forecasts of the climate-driven changes in butterfly species richness. The expected future conditions imply widespread losses of regional richness, even under the less severe scenarios. Since the negative impact of warming is expected to be extensive, long term conservation plans should concentrate in the present protected areas of highest richness as these are most likely to represent the last refuges for mountain species.


Introduction
Global warming impacts the geographic range of organisms and their regional richness patterns (Díaz et al. 2019;Cardoso et al. 2020;Wilson and Fox 2021), though to different extent across taxonomic groups and geographic areas (e.g., Aragón et al. 2010;Devictor et al. 2012;Outhwaite et al. 2020). One of the ways to incorporate expected climate change effects to long term conservation plans is to geographically forecast the diversity and composition of species assemblages under future scenarios (e.g. Hannah et al. 2007;Schmitz et al. 2015;Fordham et al. 2016), hence forecasting methodologies represent an active research topic in the study of climate change consequences on biodiversity (e.g. Drew et al. 2011;Mouquet et al. 2015;Petchey et al. 2015;Urban et al. 2016;Lovejoy and Hannah 2019).
Climate-based expectations for future decades are negative for an important number of European butterflies (Lepidoptera, superfamily Papilionoidea) (Settele et al. 2008;Warren et al. 2021). Butterflies are good indicators of environmental variation including that of climatic nature, so any mismatches in their interactions with their resources may rapidly translate to shifts in their spatial distribution (Wilson and Maclean 2011;Schweiger et al. 2012). In fact, 1 3 the negative expectations at the mid-to broad scale in the area are consistent with the negative population trends found for most species in the western Mediterranean Basin (Stefanescu et al. 2011;Melero et al. 2016;Colom et al. 2020): a matter of concern, since the Mediterranean peninsulas host an important fraction of the European butterfly diversity (including the Iberian Peninsula, a biodiversity hotspot: Myers et al. 2000). How the described demographic trends would eventually translate to species richness and how such changes would distribute across this specific region, remains speculative. The mid-term (e.g., decades) responses of the geographic ranges to the expected climate warming has been modelled for a few Iberian species (e.g., at grid sizes of 10 × 10 km: Romo et al. 2014;García-Gila 2019) to suggest dominantly negative expectations, but an estimate of the magnitude and geographic distribution of the net expected change in species numbers (richness) is not yet available.
Our purpose here is to estimate the amount and the spatial distribution of the potential change in butterfly species richness under future climate scenarios in the Iberian Peninsula under a macroecological approach (that is, modelling species richness, i.e., species counts per area unit, as a variable: Whittaker et al. 2001;Pereira et al. 2013). Compared to stacking individual species models (S-SDMs) this is less computationally intensive far less sensitive to the problem of modelling the localized endemic species (e.g., Pineda and Lobo, 2012; for further argumentation on the relative merits of these see e.g., Guisan and Rahbek 2011, Calabrese et al. 2014, Distler et al. 2015, Biber et al. 2019). However, forecasting distribution data to future climate scenarios requires a trained statistical model of contemporary richness whose performance is limited by a chain of methodological constraints and decisions including not only the statistical approach (e.g. Botkin et al. 2007;Guisan and Rahbek 2011;Calabrese et al. 2014;Mouquet et al. 2015) but also data quality (Hortal et al. 2007), their environmental representativeness (e.g. Hortal et al. 2008a;Sánchez-Fernández et al. 2021) and the bearing of non-climatic factors whose effects could be mistaken with, or interact with the climate predictors (Guisan and Thuiller 2005;Acevedo et al. 2017).
Among the items mentioned above, we shall concentrate in two that are merely relative to the basic information needed, i.e., the quality of the data (since this is not uniform across the study area, which constrains the geographic resolution and determines the selection of area units, as detailed in the methods section), and the contribution of non-climatic variables. The last perspective is justified by the fact that a realistic forecast of future richness should rely on all its main drivers, including those that do not directly reflect climate conditions but that may interact or be correlated with them (Luoto and Heikkinen 2008;Titeux et al. 2009;Nieto-Sánchez et al. 2015;Herrando et al. 2019). Some of those variables are temporally dynamic as the result of human activity and are difficult to use in forecasting due to the temporal extent of the available data (but cf.: Bouwman et al. 2006;Reidsma et al. 2006), while other are 'static' (i.e., virtually constant along a period of decades, such as altitude or the geological substrate). Forecasting species richness using only climatic predictors may overestimate real changes, while using solely static non-climatic variables should render forecast patterns of constant richness. Thus, results from a combined model (climate plus static variables) might generate intermediate -and perhaps more realistic-alternatives to the pure-climate based ones, largely depending on the magnitude of the interactions between the two types of variables. Moreover, the residuals or non-explained variability derived from these models may show a geographically distributed pattern which probably represents unknown predictors not accounted for by the model (e.g.: Birks 1996; Hortal et al. 2004Hortal et al. , 2008b. On this basis the spatial information may be incorporated as a part of the 'static' factors to, in the absence of other non-climatic predictors, provide a reasonable (though hypothetical) contrast of 'minimum change' to purely climate-based forecasts, helping to identify the potential range of expectations on richness change under the future climate scenarios.
Considering the sampling intensity carried out historically in the Iberian Peninsula, this study aims to provide different estimates of the range of expected change of butterfly species richness from projections to four future climate scenarios, based in models relying in different combinations of variables, both climatic and static. The intensity and the geographic distribution of the expected changes of butterfly species richness will be discussed in terms of the explanatory ability of the models fitted and their congruence with recent evidence from butterfly demographic trends, with emphasis in the sources of uncertainty of the forecasts that are of practical relevance.

The data
The study area was the Iberian Peninsula (the continental territories of Spain and Portugal) and the taxon the superfamily Papilionoidea, represented in this area by six families and 238 species. The butterfly data come from García- , updated to cover the period 1900-2019 (mean year = 1990, SD 22, ca. 400,000 records) from the literature, collections, field observations and available Internet repositories.

Grid size selection and cell completeness
Sampling effort was evaluated at five grid resolutions (MGRS squares with sides of 10, 30, 50, 100 and 200 km), based on incidence data where the rows were the unique combinations of species x coordinates x date x observer. A rational function was fit to accumulative (randomized) relationship between species numbers and records (with the package KnowBR under RWizard: Guisande et al. 2014;Lobo et al. 2018;Guisande and Lobo 2018). Sampling effort was measured as completeness (estimated as S est /S obs , where S obs = number of species recorded and S est = richness value at the curve's asymptote). The frequency of high completeness scores was poor at the highest resolution tested (10 km squares) where less than 10% of the cells reached 90% completeness (Fig. S1). The coverage of well-surveyed cells increased with cell size, but only a very large cell size (100 × 100 km or more) ensured full coverage; the intermediate 50 × 50 km grid was selected as a reasonable balance between precision and resolution. At this grid size 17 cells were excluded for containing either less than 25 species, 100 records or a records/species ratio lower than 3.0. The remaining 224 workable cells had completeness values of 87.6% on average (SD 10.8), higher than 50% in more than 90% of them and above 80% in 175 cells (though exceeding 95% in only 69 units). The geographic distribution of completeness, observed species numbers, and estimated richness at this grid size are presented in Fig. S2.
We initially used the 80% threshold for completeness to select the 'reliable' cells for modelling (comparable to e.g., Pelayo-Villamil et al. 2018;Lobo et al. 2018). However, the area units selected for modelling should cover the range of the main environmental gradients, while collecting biases are frequent in chorological data sets (typically, more records from more species-rich areas: Hortal et al. 2008a;Sánchez-Fernández et al. 2021). This problem (which is likely to affect our data: Romo et al., 2006) would impose differential collecting intensity along at least one environmental gradient to result in a part of that gradient being excluded from the training data set. Alternatives to solve this may include modelling on the whole data set while controlling for completeness, either weighting values proportionally to sampling intensity, or including sampling intensity as a covariate. Thus, we tested the three alternative procedures ('good' cells, case weighting, and entering completeness as a covariate) in combination with the remaining options described below. Contemporary (1970Contemporary ( -2000 climate data were derived from the data from WorldClim version 2.1 (30'' resolution: Fick and Hijmans 2017). Nearly all the WorldClim variables (BIO1-BIO18) are reciprocally correlated in the study area; after a preliminary selection we retained the following nine: mean yearly temperature (BIO1), maximum temperature, average and absolute maximum (both from BIO5), minimum temperature (mean and absolute minimum, both from BIO6), temperature annual range (BIO7), mean annual precipitation (BIO12), mean annual range of precipitation (the difference between BIO13 and BIO14, precipitations of the wettest and the driest month) and summer precipitation (the monthly precipitation in the warmest quarter, BIO18). The static variables included altitude (derived from ETOPO2v2: National Geophysical Data Center 2006), lithology (from the 1:200,000 geologic chart by IGN, 2018), plus elementary distance or position metrics derived from the data coordinates: mean, minimum and maximum altitude and the altitudinal range; the percent coverage of acid rock, basic rock and clay / sediments, the lithological diversity (variance from the three former variables), the balance of acid / basic rock (on the scale 1 = 100% acid to 3 = 100% basic). Finally, Latitude and Longitude (the cell centroids MGRS coordinates, transformed to decimal degrees), the shortest distance to the coast and the surface of the cell occupied by land (Also see Supplementary file S1).

Sampling effort and dominant environmental gradients
The representativeness of the grid cells in terms of the main environmental gradients was assessed from the first two factors from independent PCA analyses applied to the full set of cell values for the variables describing contemporary climate, physiography, and geographic position.
The factors values were classified into 10 equal intervals and the cell scores were plotted against the cell's completeness for visual inspection.

Shared explanation in the test variables
Variance Partitioning (e.g., Peres-Neto et al. 2006) was used to identify the a priori proportion of shared effects of the climate, physiography and position variables on richness variation as well as to explore some of the modelling results was calculated with SAM (Rangel et al. 2010).

Model fitting
The explanatory variables were standardized. Generalized Linear Models (GLM) were applied assuming a Poisson distribution of richness and a log link function with package 'car' (Fox and Weisberg 2019) under R vers. 3.6.1 (www: cran.r-project.org). The variables were tested in descending order of their bivariate correlation with richness (Table S1). Variable selection was done manually on a forward-backward stepwise sequence. Deviance and autocorrelation were checked at each step. Autocorrelation was measured via the variance inflation factor (VIF: Fox and Monette 1992) with a threshold VIF > 3. Backward selection was performed after any new variable was added, and at that point any redundant variables or those inducing VIF values above the threshold were eliminated. To account for curvilinear relations, the quadratic terms of the variables were tested and added if the linear monomial of the same variable remained in the model. Model fit was measured by adjusted deviance reduction (D 2 adj : Guisan and Zimmermann 2000). Validation was done through the Mean Square Error (MSE) from the 'leaveone-out' cross-validation procedure (package 'boot': Canty and Ripley 2020) with 45 cells of completeness higher than 97% as the validation set.
After each model was fit, the spatial structure remaining in the residuals was measured with Global Moran's I (Diniz-Filho et al. 2003) at distance lags of 50 km (with package 'ape': Paradis and Schliep 2018) as well as performing a Trend Surface Analysis (TSA: Legendre and Legendre, 2012). Any other calculations were done with Statistica (StatSoft Inc. 2004).
Alternative models 18 models were tested, differing in the choices: (i) climate-only models vs. climate plus non-climatic variables models; (ii) how sampling effort was incorporated into the models, and (iii) alternative selection of climate variables. The combinations defining each model were abbreviated with three capital letters in the order (i-ii-iii). More in detail: (i) Climate vs. climate + static variables, A/C/K where: A = any of the available variables was selected, C = only the climate variables were tested and K = like A but with the highest possible proportion of the variance attributed climate. To do this, all the static variables were regressed onto the climate variables entered in each model and their residuals taken. These residuals (representing information not included in the C models and linearly independent from the climatic variables) were tested sequentially and added to the former C models if their effects were significant.
(ii) Incorporating cell completeness as a measure of realized sampling effort in either one of three ways (S/W/C): S = only the cells where completeness > 79.9% were analyzed (n = 175). W = the estimated sampling effort was used to weight the cases. Since both the slope of the accumulation curves and completeness equally describe the accuracy of S est and these were inter-correlated (r = − 0.89, P < 0.001) we calculated their geometric average [((1-slope)·(comple teness))^0 .5 ] and re-scaled the result to the range 0.1-0.9; these values were incorporated as weights to the cases in weighted regression. And C = Completeness was entered as a co-variable in the models. (iii) Including or not summer precipitation (a/b). a = models were fitted with no restrictions on the climatic variables, or b = summer precipitation was omitted. This derived from the finding that in some cells under future climate scenarios (namely 2070-2.6) summer precipitation might increase despite of a decrease of the total rainfall (Fig. S3). The sense of this change is unprecedented and might unveil summer precipitation as an equivocal predictor of future richness, because in the contemporary estimates the total precipitation and summer rainfall are correlated all over the area.

Richness change in future scenarios
The values of the future climate variables were transformed using the means and standard deviations of the corresponding contemporary data and then used to forecast from the models selected. Richness changes were measured as the absolute differences (future-present) between the cell forecasts and the contemporary values, the latter being represented both by the estimated ones (S for − S est ) and by the model predictions (S for − S pred ). The 'certainty' of the richness shifts per cell was measured as a/b for the difference (S for − S est ) and as a/(bc) for (S for − S pred ), where a = the variation from model averaging (s.d. within one scenario), b = the error associated to present estimates of richness (measured as S obs /S pred , which is completeness), and c = the raw cell residuals from the regressions.
The set of cell forecasts was compared by one-way ANO-VAs with the model types and climate scenarios as factors. The same data were submitted to a PCA analysis from which the two first factors facilitated a graphic comparison of the relationships between the model types and the cell results.

Sampling effort, cell completeness and dominant environmental gradients
The variance accounted for by the first two components from each of the three subsets of variables (PCA) was 98.5% (position), 65.3% (physiography) and 75.2% (climate). Within these limits at least one cell with completeness > 90% was present in each section of the ranges of any component (Fig. S4, Tables S2 and S3). Richness was significantly correlated with temperatute (mean temperature, r = 0.78), maximum altitude (r = 0.74), summer precipitation (r = 0.71) and less so with the lithological nature of the substrate. Completeness was significantly correlated (r = 0.20 to 0.40) with species richness and with several predictor variables including those already mentioned (details in Table S1). Variance partitioning showed a good explanatory capacity by a combination of the candidate variables but a low 'pure' contribution by each subset (less than 1% in the case of climate variables: Tables S3, S4 and S5).

Models selected
The models combining climatic and static variables offered the best fit in terms of Deviance (D 2 ), followed by K models (climatic variables with residual effects attributed to the static variables); climate-only models had the lowest performance (full details in Supplementary Tables S6, S7 and S8). In the combined models, the shared proportion of the explained variance was meaningful (30% to 47%) and only in two of them (Table S7: ACa, ASb) the 'pure' climatic part of the explanation exceeded 40%. At least one measurement of temperature and one of precipitation tended to enter in every model (though temperature was replaced by elevation in ASa and AWa). Omission of summer precipitation resulted in the selection of the minimum temperature rather than e.g., the average precipitation. A comparison of all the models (through PCA analysis on the whole set of forecasted richness values: Fig. S5) suggested two slight general trends: 'climateonly' to 'all variables' models, and 'with-to without summer precipitation' ones.
The spatial structure of model residuals (Moran's I, TSA) tended to be small (less so in the climate-only models). Between 64 and 94% of the initial values (I = 0.231, P < 0.001) was accounted for in the results (except for model CSb: Table S6). Interestingly, the lowest remaining spatial effects occurred when completeness was included as a covariate.

Richness change in future scenarios
The average cell expectations for species richness change are shown in Fig. 1 (after discarding models with completeness as covariate, given the correlation between completeness and other explanatory variables). Figure 2 shows that the magnitude of the estimated species richness change was proportional to the intensity of change in the climatic scenarios selected (i.e., from most moderate in 2070_2.6 to larger in 2070_8.5) and summarizes the direction of the changes imposed by the model building choices (further details are provided in Supplementary file 1: Tables S9, S10, S11 and S12 and Fig. S7; and detailed results by cell and scenario in Supplementary file 2).
Overall, species richness shifts were highest in the less thoroughly surveyed cells (r = − 0.346, P < 0.001) and for cells with the highest raw model residuals (r = 0.519, P < 0.001), as in fact the raw residuals and completeness were negatively correlated (r = − 0.327, P < 0.001) (Figs. S5, S8). In other words, the most favorable forecasts tended to correspond to the less thoroughly surveyed cells and are associated to comparatively poor model fit for those cells.

Data quality, geographic and environmental coverage
Despite a long history of Iberian butterfly faunistic studies (see summary in García- ) our data did not provide a dense coverage of well-surveyed cells at a reasonably high resolution such as 10 km cells. The loss of resolution was compensated by the high coverage of the environmental gradients by the well-surveyed cells, suggesting that extrapolation from heavily inaccurate species richness values was avoided (as emphasized by e.g., Hortal et al. 2007 andRocchini et al. 2011). Our relatively poor results on completeness show that Iberian butterflies faunistics still have a way to go. Moreover, the historical accumulated butterfly data probably reflect a degree of collecting bias (Romo et al. 2006) as shown from other European countries (Dennis and Thomas 2000). Considering the rapidly evolving environment and the limited resources devoted to the collection of primary data, a planned recording scheme is needed (as already pointed out by e.g., Sastre andLobo 2009 andSánchez-Fernández et al. 2022).

Incorporating sampling intensity to the modelling process
Judging from the best-fit statistics, trusting the bestsurveyed cells was the best alternative to deal with completeness, at the cost of shortening sample size. Adopting sampling effort as a covariate was less efficient and, in our specific case, led to less reliable predictions because of the correlation between completeness and other predictor variables. Weighted regression offered the poorest D 2 scores. We might simply conclude that relying on the best-known cells is the evident choice (e.g.: Pelayo-Villamil et al. 2018;Ronquillo et al. 2020), but a cautionary note is pertinent here: if the recording patterns behind our data base are biased towards the most speciesrich areas in coincidence with some of the environmental gradients ultimately explaining butterfly richness, best-fit models may be not completely realistic (further discussion below) even when there is a priori evidence for full coverage of the main environmental gradients. Fig. 1 Mean expected difference in richness (future-present) by reference to present estimated values (from accumulation curves, left column) and to present predicted values (predicted by the models, right column). The values shown in the legend are the limits of the ranges used for colour labelling. Certainty (see text) is represented by black dots. All values averaged from the 12 models where Completeness was not used as a covariable (details in Tables S4, S5 and S6). The cells marked with white circles are those without sufficient information to calculate estimated richness; this was replaced by their know species numbers (Sobs) for graphic presentation only, their values may be not reliable

Good correlations and false cognates
Regarding the performance of the alternative subsets of variables to model contemporary species richness, the combinations of climatic with the non-climatic variables rendered the best model statistics. However, the small 'pure' contribution of each subset of variables indicates that identifying the causality of richness from our data set is not straightforward (Lobo et al. 2002 and discussion therein;Legendre and Legendre 2012). The risks of forecasting temporal changes from spatial variation (Kerr et al. 2011) increase if the future spatial and temporal relations between the variables vary beyond the range tested (Araújo and Rahbek 2006;Dormann et al. 2012). For instance, replacement of summer precipitation by alternative candidate variables increased model fit, despite the initially highest correlation between that variable and species richness. Including or not this variable imposed the most significant differences among model predictions: omitting it increased model fit and resulted in more negative expected richness shifts. Most of the Iberian Peninsula is dominated by a Mediterranean climate featured by spring and autumn rain periods combined with summer drought (Aschmann 1973;Belda et al. 2014). An inverted ratio of summer to yearly precipitation is expected, but it is unlikely that an increased summer rainfall combined with an overall yearly water shortage would enhance the potential for richness of Mediterranean-adapted organisms. However, butterfly richness is positively correlated to summer precipitation in the area (perhaps via its correlation with productivity and ultimately with actual evapotranspiration: Aragón et al. 2019), which makes this variable a misleading predictor of future richness.

The intensity and geographic distribution of expected richness changes
The species richness changes expected under the future climate conditions tested were negative overall, more severely under the most extreme climate scenarios (consistently with the increasing temperature and decreasing precipitation expectations in the Mediterranean: Giorgi and Lionello 2008;Hertig and Jacobeit 2008). This makes sense because in this area butterfly richness is positively associated with altitude and precipitation and negatively with temperature (Romo et al. 2006; and our results). The shifts (future-present) are best understood if the comparisons with contemporary richness represented by model-predicted values are considered first (Fig. 1, right column; compare with Fig. S2).
The results indicate generalized losses (94% to 98% negative cell values, average losses above 11 species per cell), with only a small area in the central Pyrenean mountains (1 to 3 cells depending on the emission scenario) showing neutral or slightly positive changes. The heaviest richness  Table S5) with 95% confidence limits. All the differences between means within each of the four factors are significant (p < 0.01, post-hoc Fisher tests) except for the pairs connected by dashed lines and labelled 'n.s.' loss should concentrate in the inland eastern areas in the shortest term, proceeding to profound negative changes in the central-eastern sectors or virtually across the peninsula (depending on the emission scenario), with comparatively moderated effects (though still negative) along a 50-100 km stripe along the Atlantic coasts (N, NW and SW). The expected variation is slightly more benign by reference to the estimated values (S est ), though still indicating a mean loss of 10-28 species per cell and negative values in 71% to 87% of the area. Under this comparison (Fig. 1, left column) the dominant negative values are again combined with less severe losses and even some positive values; these spread along the northern and western coasts but also occur in scattered inland cells, more frequently in the SE half of the peninsula under the less intense emission scenarios (2050-2.6, 2070-2.6). A few of the 'positive inland spots' persist under the most adverse 2050-8.5 and 2070-8.5 conditions.

Relevance of the standards for contemporary richness in incompletely surveyed areas
The finding that the forecasts included some positive cell values -opposite to the negative average trend-is remarkable. The forecasts differed markedly depending on whether contemporary species richness was represented by the model predictions or by the values derived from accumulation curves. This must reflect environmental variation not accounted for by the models and is likely to result from overpredicted contemporary richness in species-poor cells (also see Hortal and Lobo 2011;Calabrese et al. 2014;Biber et al. 2019). This is consistent with our results (completeness is positively correlated with richness and negatively with the model residuals) and explains why forecasts for areas of modest richness often led to positive richness shifts only by reference to the fitted values. Two complementary interpretations are possible. The first is that some sources of environmental variation were not represented by our predictor variables and are correlated with species richness and sampling intensity. The second is that, with descriptors for human-driven effects missing among our variables, the inflated richness values for some cells do in fact represent the contemporary effects of human activity, hence an extinction debt (Malanson 2008;Figueiredo et al. 2019).
An optimistic interpretation of these results (i.e., the expectation for real local increases of richness) is not straightforward. Our forecasts must be interpreted in macroecological terms, and from this point they fit the expectations in terms of the relationship between richness and the water-energy balance or productivity Whittaker et al. 2007). In this context, positive cell forecasts represent an enhanced potential for richness from species similar to those from which the richness patterns were modelled; this might be realized, or not. Most of the 'positive' cells are coastal ones or stay outside the mountain ranges, while contemporary richness maxima occur in the mountain ranges. Range spreading of the mountain-adapted species is unlikely in a scenario dominated by raising temperature, and northward immigration from North Africa is largely prevented by the Mediterranean Sea. So, richness gains would most likely result from the spreading ranges of one part of the current moderately widespread and resilient species (cf. Melero et al. 2016); in other words, via changes of the geographic turnover of species composition and perhaps implying a degree of faunal homogenization (as recently reported from French bird communities: Rigal et al. 2022).

Potential role of the static variables
The forecasts derived from our climate plus static variables models deserve a comment. We are aware that climate variation represents only one part of the changes expected to impact biodiversity along the next decades, but that no dynamic variables other than the climatic ones were explicitly modelled. Instead, we kept the constant effects of non-climatic and positional variables in the forecasts to attenuate the warming effects (for instance, recall that spatial position by itself alone provided a reasonable explanation of richness). Although the results were in the direction predicted, the intensity was small or negligible. Our results may partly cover the potential buffering effect of topographic heterogeneity to warming conditions for butterflies (Oliver et al. 2014;Suggitt et al. 2018), because topographic variation strongly co-varies with altitude in the Iberian Peninsula. There is room for suspecting that if any variables not dealt with here have a bearing on richness, are correlated with physiography or have a spatial structure and are likely to change in the near future, richness shifts should be stronger than indicated by our results and most likely negative. Recent evidence from butterflies in Spain demonstrates changes in demography, species composition and altitudinal ranges running parallel with warming, land management changes, or both (Wilson et al. 2007;Herrando et al. 2016;Ubach et al. 2019;Warren et al. 2021). Land use in the study area has undergone profound changes along the last century in parallel to warming, both probably affecting synergistically butterfly demography (as stated by Stefanescu et al. 2011) to complicate the identification of the relative impacts of both types of change. In fact, on the wide geographic scope, impacts to Mediterranean biodiversity are expected not just from warming but from a combination of threats of which human land use changes outcome the relevance of climatic variation (Sala et al. 2000).

Iberian trends of butterfly richness and the chances for preservation
To conclude, our results on butterfly species richness are broadly consistent with the range contractions predicted for an important number of species in SW Europe (Settele et al. 2008;Romo et al. 2014) which would be the ultimate result of the widespread specific negative local demographic trends already documented: the climatic perspectives for the Iberian butterfly fauna within the next 50 years are negative through most of the area. Moreover, compensation of the negative impact of warming by non-climate variables may be limited unless these come from variables completely uncorrelated with the ones tested here.
The present network of protected areas appears to be efficient for hosting virtually the full set of Iberian butterfly species including the endemic ones (Romo et al. 2007;Rosso et al. 2018). However, as a rule, protected areas have been reasonably surveyed for butterfly diversity and tend to be species-rich: with few exceptions, these are just the areas where climate change will impact more severely after our results. Although no precise indications can be extracted from our coarse resolution maps, special attention should be paid to improve the potential for butterfly richness of the highest diversity areas within the present network of preserved sites, as some of these may ultimately act as the last refuges for some Iberian populations of mountain butterflies. Concern on the limited efficiency of the available corridors between the Iberian protected areas Lobo 2018, 2021) is more than justified for the butterfly fauna because of the unfeasibility of interconnecting mountain areas. Efforts to re-assess the contemporary drivers of butterfly richness at a finer scale and a short-term planned recording scheme to complement the existing distribution data are needed before the expected impacts of change on a finer geographic resolution prove necessary.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.