Potential distributions of invasive vertebrates in the Iberian Peninsula under projected changes in climate extreme events

Invasive alien species (IAS) can cause profound impacts on ecosystem function and diversity, human health, well‐being and livelihoods. Climate change is an important driver of biological invasions, so it is critical to develop models and climate‐driven scenarios of IAS range shifts to establish preventive measures. In this study, we analyse how projected changes in the frequency and magnitude of climate extreme events could affect the spread of the six most widely distributed invasive vertebrate species in the Iberian Peninsula.


| INTRODUC TI ON
Invasive alien species are among the five direct drivers of environmental change with the largest relative impacts on biodiversity and ecosystem services (Brondizio et al., 2019;Early et al., 2016). They can also affect human health through the introduction and spread of new pathogens, and they can jeopardize human well-being and their livelihoods by causing infrastructure damage, food destruction and reduced economic growth due to both decreased performance of human activities and increased expenses to reverse or mitigate the losses they cause (Pyšek et al., 2020). Estimates indicate that invasive species have cost the European Union (EU) at least 20 billion Euro per year over the last decades (Kettunen et al., 2008). Actual costs are probably higher, as many countries have begun to account for such costs only recently (Jardine & Sanchirico, 2018).
Although the number of invasive alien species has been rising over time worldwide (Seebens et al., 2017), Europe, in particular, has become a hotspot (Pyšek et al., 2020), with circa 14,000 alien species and roughly 10% of them being considered invasive according to the European Information Network of Exotic Species (EASIN).
This increasing trend is likely to continue in Europe in the coming decades, as rates of biological invasions are predicted to remain high in economically developed countries as climate change unfolds, and tourism and plant and pet trade intensify (Early et al., 2016;Ribeiro et al., 2019;Seebens et al., 2021).
Climate change is a key factor that might underlie future movement of alien species (e.g. Bellard et al., 2013;Peterson et al., 2008).
On the one hand, some previously harmless alien species are becoming invasive as climate change boosts their competitive advantages with regard to native species (Pyšek et al., 2020). On the other hand, several species are also shifting their natural ranges without human intervention in response to climate change and thus expanding into regions in which they previously could not survive and reproduce (Peterson et al., 2008;Walther et al., 2009).
Nevertheless, there is a lack of consensus on climate change's net effect; while some modelling studies indicate that climate change increases the area invaded by invasive species, others conclude that it limits their distributions (Bellard et al., 2018;Gallardo et al., 2017).  . Therefore, there is an urgent need to better comprehend how climate change drives invasive species distributions (Bellard et al., 2018), so that preventive measures can be established accordingly. For instance, Orlikowska et al. (2016) suggested that lack of modelling studies limits insight enabling adequate management of potential impacts of invasive species in protected areas. Although there is no shortage of climate change modelling studies in general, there is a generalized lack of empirical evidence on the reliability of their predictions (cf. Barbet-Massin et al., 2018). Many experts are also highlighting the need to develop models and future scenarios of biological invasions at the continental scale for horizon scanning (e.g. Essl et al., 2019). However, it is important to bear in mind that forecasted effects (direction and magnitude) of climate change on invasive species distributions at large spatial scales often differ from those forecasted at local or regional scales, where most management actions are undertaken (Bellard et al., 2018).
The changing climate is leading to changes in the frequency, intensity, spatial extent, duration and timing of climate extremes (Fonseca et al., 2016;Vicente-Serrano et al., 2014), which in turn are driving abrupt ecological changes resulting in range shifts in species from a variety of taxa (Ummenhoffer & Meehl, 2017). In particular, extreme weather or climatic events such as heat waves, droughts, heavy precipitation, floods or tsunamis may facilitate the spread and establishment of invaders by (a) allowing them to reach further distances (Čuda et al., 2017), even new regions (Carlton et al., 2017), (b) opening new niches that they can quickly fill (Straub et al., 2019) and (c) creating stressful conditions leading to reduced competition and predation that diminishes the biotic resistance of communities (Vetter et al., 2020). (See review by Diez et al. (2012) for a deeper understanding of these invasion mechanisms linked to climate extremes.) Invasive species are usually generalists and tolerant to a wide range of conditions, and they tend to spread fast and thus have the capacity to rapidly colonize disturbed habitats; these traits may allow them to cope and thrive during and after extreme climate conditions, thus out-competing native species (Diez et al., 2012). By contrast, certain invasive species that spread and replace native species during periods of moderate climate may collapse under extreme climate events, which might compromise the ecological resilience of the invaded ecosystems (Alpert, 2006). In general, whether a certain extreme climate event that increases environmental stress levels favours either native or invasive species is taxon-and system-specific (Diez et al., 2012). These processes might be particularly relevant in the Mediterranean region, wherein extreme heat waves are projected to become more frequent, more intense and longer-lasting, leading to increased frequency and intensity of both droughts and heavy precipitation events (e.g. Prasad Dasari et al., 2014); these risk for the species examined although increases in favourability should be expected locally.

K E Y W O R D S
biological invasions, climate change, conservation planning, distribution shifts, extreme climate events, invasive alien species, species distribution models phenomena, in turn, are likely to augment the length and severity of the fire season and the risk of river floods, respectively (EEA, 2017).
Species distributions modelling studies frequently overlook the effects of climate extremes on spatiotemporal invasion dynamics.
To address this issue, we analyse how projected changes in extreme climate events might affect the spread of various invasive vertebrate species in the Iberian Peninsula. We focused on six well-established, widely distributed and listed as invasive species of concern in both Spain and Portugal. These include four bird species: Red avadavat Amandava amandava Linnaeus Common waxbill Estrilda astrild Linnaeus Monk parakeet Myiopsitta monachus Boddaert, and Roseringed parakeet Psittacula krameri Scopoli; one mammal species, the American mink Neovison vison Schreber; and one reptile species, the Pond slider Trachemys scripta Schoepff. To forecast potential favourable areas for invasion, we built ensembles of species distribution models (SDMs) (Araújo & New, 2007) using, as predictor variables, extreme climate events metrics projected onto the future based on different climate change scenarios (see section 2.2 in Methods for details). SDMs characterize species' current potential distributions, that is the climate space where these invaders have already been found, and project changes in the distribution of climate space according to future climate change scenarios .
This process enables early detection of areas with high risk of invasion, thus allowing a rapid response where it is most necessary and effective. Given the projected changes in climate extremes for the Mediterranean region, we expect a reduction in the favourable area of invasive species exhibiting low tolerance to extreme warm temperatures and that require a markedly wet season to survive and grow.

| Study area and species distribution data
The Iberian Peninsula is located in the southernmost western tip of Europe and includes the mainland territories of Portugal and Spain, covering an area of ca. 600,000 km 2 . The climate and topography are heterogeneous and the region has a marked peninsular character, as its connection to the European continent is relatively narrow (about two-fifths of the northern boundary) and is crossed by the Pyrenean mountain range. This limits both biotic and abiotic interactions with adjacent territories, making the Iberian Peninsula a discrete unit adequate for biogeographic analysis.
Species occurrence records were gathered from the Spanish and Portuguese distribution atlases of birds (Equipa Atlas, 2008;Martí & Del Moral, 2003), mammals (Bencatel et al., 2019;Palomo et al., 2007) and reptiles (Loureiro et al., 2008;Pleguezuelos et al., 2002). The spatial resolution of grid cells was 10 × 10 km 2 . These records were complemented with additional information from reference books and peer-reviewed scientific literature (e.g. Ascensão et al., 2021), as well as interviews with wildlife managers and experts for each taxonomic group (see Appendix S1 for details). We also added records from online databases, such as the Global Biodiversity Information Facility (GBIF.org, 2020) and the European Information Network of Exotic Species (EASIN, https://easin.jrc.ec.europa.eu/easin). From these records, we excluded the ones whose occurrence status was reported as "absence," whose abundance (or individual count) was reported as zero, or those whose reported coordinate uncertainty was larger than the spatial resolution of our analysis grid.
Additionally, we used the "scrubr" R package (Chamberlain, 2020)  Although many occurrence records were at a finer spatial resolution, for the sake of standardization we used the same 10 km × 10 km (i.e. 100 km 2 ) resolution for the entire analysis.
Regardless of the total number of observations per cell, each occupied cell counted as only one presence. This is analogous to thinning occurrence data at 10-km intervals, and it reduces the influence of sampling bias due to some especially accessible localities having disproportionally large survey efforts (e.g. Pereira et al., 2020).
Any spatial autocorrelation left in the data, that is 100-km 2 cells having increased probability of presence because they are close to other 100-km 2 cells with presence, is a natural result of dispersal processes and should not be artificially removed (Legendre & Legendre, 2012). Model residuals were obtained for each model by subtracting to the observed values the predicted values obtained with the "fitted" R function (using type = "link" to get them in the scale of the predictors). The spatial autocorrelation of these residuals was then calculated through Moran's I statistic with the "moran. mc" function of the "spdep" R package (Bivand et al., 2013) using 999 random permutations. We further dealt with the influence of these spatially contagious processes by factoring in the spatial trend in each species' current distribution (see section 2.3.2) and by selecting models based on spatial block rather than random crossvalidation (see section 2.3.1).

| Climatic variables
To characterize extreme conditions, we computed nine variables based on temperature, precipitation and humidity (Table 1) Table S3.1.
We used RCMs instead of global climate models (GCMs) because the latter have an overly coarse resolution compared to the spatial resolution of our species distribution data. As RCMs downscale climate fields from GCMs, they provide information at fine (meso or micro) scales that are more accurate for studies of regional phenomena and for application to climate impact assessments (Giorgi, 2019;Jacob et al., 2020). GCMs, in contrast, have been found to poorly represent actual climatic conditions in a significant part of our study area, a problem that is exacerbated when models are used for future projections of species distributions (Bedia et al., 2013).  (Figure 1), although the final predictions for these cells were interpolated based on squared inverse distance, using the "idw" function of the "gstat" R package (Pebesma, 2004). Spearman's Rank Order Correlation analyses were performed and highly correlated variables (Spearman's rho coefficient > 0.75) were removed to avoid multicollinearity (Dormann et al., 2013).
Although previous studies have highlighted the importance of using both the native and invaded ranges when predicting the ex-  (Muñoz & Real, 2006). Models may even show higher predictive accuracy when using invasive data only (Barbet-Massin et al., 2018). While global occurrence data might allow a better appraisal of the full potential invasive range, from a management point of view it is more relevant to predict the areas that are most likely to be invaded next (Barbet-Massin et al., 2018).

F I G U R E 1
Occurrence data of the modelled species in the Iberian Peninsula (SW Europe). Presences (in black) are on Universal Transverse Mercator (UTM) grid cells of 10 × 10 km 2 . Grid cells at UTM zone boundaries (white cells) were not included in the analyses for having nonstandard sizes

| Species distribution modelling
We followed best-practice standards in species distribution modelling regarding guidelines for response and predictor variables (see previous section), and for model building and evaluation Sillero & Barbosa, 2020). We provide as Supporting Information (Appendix S2) a description of the modelling steps taken to develop our SDMs following the ODMAP (Overview, Data, Model, Assessment and Prediction) protocol (Zurell et al., 2020).
We also provide the R code used in our analysis, together with most of the data (except where redistribution was not allowed) in a public GitHub repository (https://github.com/AMBar bosa/extre me-invaders).

| Model building and validation
We used an ensemble approach with four SDM methods of varying complexity: Generalized Linear Models (GLM), Generalized Additive Models (GAM) with smoothing splines, Random Forests (RF) and Bayesian Additive Regression Trees (BART). All these methods produce presence probability, which is important for a commensurable combination of predictions from different methods and scenarios (see below). GLM, GAM and RF have been used widely in SDM (e.g. Peterson et al., 2011), while BART is a newer method producing highly promising results (Carlson, 2020; see also our results below).
With random forests, it is customary to compensate for imbalanced datasets, either by attributing equal weights to presences and absences, or by using the same number of presences and absences in each decision tree. However, we needed to use the same naturally imbalanced data (i.e. reflecting the actually observed presence/ absence ratios) across species and modelling methods, so that the resulting presence probabilities could be formally compared and combined (see section 2.3.2). To ensure that using balanced or imbalanced datasets had no effect on our main results and conclusions, we also computed the random forest models with balanced data, by setting the sample size to be the same for both presences and absences in each species. To avoid wasting data with this subsampling for balance, we increased the number of trees to 1,000, to ensure that all data points were used in at least some of the trees.
The data were initially partitioned into 10-fold cross-validation sets: the grid cells of the study area were split into 10 groups, or folds, using R package "blockcv" version 2.1.1 (Valavi et al., 2019), with a systematic selection and assignment of spatial blocks into ten 100 × 100 km 2 folds ( Figure S3.1 in Supporting Information).
Block cross-validation is more reliable than random cross-validation, which can underestimate prediction error and result in inappropriate model selection (Araújo, Whittaker, et al., 2005;Valavi et al., 2019).
Spatially separating training and test datasets allows assessing if models perform well both in nearby and in more distant areas, which is especially important if models are to be extrapolated in space or in time (Valavi et al., 2019). For each combination of modelled species, RCM and modelling algorithm, 10 different SDMs were built, in which each fold was sequentially left out of the training data and used for validation of model predictions (Fielding & Bell, 1997). This implied computing a total of 1,200 models (= 6 species × 5 RCMs × 4 modelling algorithms × 10 replicates) in the cross-validation phase.
The cross-validation performance of each species, RCM and modelling algorithm combination was assessed with R package "modeva" version 2.0 ( Barbosa et al., 2013). It has been advocated that different evaluation metrics should be used and their information integrated during model selection (Lake et al., 2020). Therefore, we used three metrics that focus on different perspectives of model that is the proportion of correctly classified presences and absences (Allouche et al., 2006), using training prevalence as the threshold value for classifying predicted presences and predicted absences; and (c) Miller's calibration slope, which assesses model reliability, that is the overall deviation of predicted probabilities from observed occurrence frequencies (Miller et al., 1991;Pearce & Ferrier, 2000).
The values of these metrics were then averaged across the 10 crossvalidation folds for each species, modelling method and RCM.
We excluded from subsequent analyses models whose average cross-validation performance failed to meet an acceptable threshold for any of the model evaluation metrics. For AUC, this performance threshold was set to 0.7, which delimits poor predictions (Rapacciuolo et al., 2012). Considering that AUC can range from 0 to 1 and TSS can range from −1 to 1, the proportional performance threshold for TSS was set to 0.4. For Miller's calibration slope, we used a performance threshold of 0.5 above or below the ideal value of 1.

| Model post-processing and projection
After cross-validation, retained models were fitted again using the complete dataset (i.e. without leaving out validation data) given the sensitivity of model projections to removal of data points . Predicted probability values were then converted to favourability, which is a function of presence probability and the presence-absence ratio of the species in the modelled sample . This removes the effect of species prevalence on the predicted values, thus making predictions commensurable and directly comparable across species and time periods (Acevedo & Real, 2012;Real et al., 2017).
Current species distribution ranges can restrict future ranges due to dispersal constraints, as ranges are likely to shift gradually from their current positions over the projected time frames (Estrada et al., 2016). Indeed, a spatial restriction of the predicted values to areas deemed reachable by the species generally reduces overprediction and improves overall model performance (Mendes et al., 2020). Therefore, for each species, we intersected final climatic favourability with the favourability resulting from the general spatial trend in the occurrence pattern, defined by a third-degree polynomial of the spatial coordinates of the presence centroids (Legendre & Legendre, 2012) (Appendix S3: Figure S3.2). The use of a spatial restriction variable can be especially important for invasive species, whose occurrence patterns are largely determined by their introduction points. These points are not always related to the environmental variables that drive the species' occurrence, but they inevitably leave a footprint on their distributions in the invaded areas (Pereira et al., 2020). Rather than converting continuous model predictions into binary predictions via the application of abrupt and often unstable thresholds (Liu et al., 2005;Nenzén & Araújo, 2011), we preserved the continuous predicted values, and we used fuzzy logic and fuzzy set theory (Zadeh, 1965) to intersect spatial and climatic favourability, using the "fuzzyOverlay" function of the "Fuzzysim" R package (Barbosa, 2015). To summarize model predictions for each species and measure their consensus, we used the "prcomp" function of R "stats" package to perform a principal components analysis (PCA) among the favourability predictions given by the different modelling algorithms and RCMs. The first axis of this PCA captures consistent spatial patterns in predicted favourability across the different models Araújo, Whittaker, et al., 2005;Marmion et al., 2009;Thuiller, 2004). As favourability values are already on the same scale for all species (Acevedo & Real, 2012), the PCA was performed without standardizing input variables. We then summarized the favourability predictions for each species by calculating a weighted mean of the predictions across models, using the loadings of the first PCA axis as a weighting factor. In this way, the final prediction for each species incorporates the degree of consensus among models, which dictates how much weight each model has in the prediction, thus avoiding disparate predictions to be blindly mixed and averaged out.
Favourability models were projected into the future, using the "predict" function of R "stats" package, to forecast favourable areas for invasion by each species under the projected changes in extreme climatic events. Consensus among the predictions of models applied to the future climate was also assessed through a PCA of the favourability values, using the same procedures described above. Finally, we measured and mapped the difference between current and future favourability for each species, and the significance of that difference using paired t-tests. We also quantified overall change (gain, loss and maintenance) in favourability between the current and future periods, using the "fuzzyRangeChange" function of the "Fuzzysim" R package (Barbosa, 2015). In addition to the weighted mean of favourability predictions, we also computed their fuzzy intersection, which corresponds to the minimum favourability across models for each species (Zadeh, 1965). This means that, whichever modelling algorithm or RCM is closer to reality, favourability will be at least this one for each species. Favourability maps were produced using packages "rgdal" (Bivand et al., 2019), "rcolorbrewer" (Neuwirth, 2014) and "cartography" (Giraud & Lambert, 2016).
Projecting models to different environments from those where they were built often entails extrapolating into novel environmental conditions, that is beyond the range of values that the models were able to analyse. In those conditions, model projections need to be taken with additional care. We assessed novelty in the projected environmental spaces using multivariate environmental similarity surfaces (Elith et al., 2010) calculated for each RCM with the "MESS" function of the "modeva" package.
BART models showed the best balance regarding discrimination, classification and reliability, and all 30 BART models met all crossvalidation thresholds for every species and RCMs (Appendix S3:

| Model predictions and projections
The first axes of the corresponding PCAs (PC1s) Figure S3.21).
The overall difference between current and future climate favourability values, averaged according to their weights in PC1s ( Figures 2 and 3), was significant for all species (paired t-tests, p < .001). Most species were predicted to experience an overall negative balance between favourability gains and losses in the study area, even if the difference was small (Figures 4 and 5). That is, the total area with favourable conditions would not increase and current unfavourable areas would mainly remain so for most studied species under the future climate. Net losses were projected to be largest for N. vison, followed by E. astrild and P. krameri ( Figure 4). However, these projected net losses do not simply reflect potential range contractions in the study area, but rather spatial shifts in favourable areas. Indeed, while some regions may become less favourable, several other regions were projected to have higher favourability for most studied species under the projected future climates ( Figure 5).

| Predicted changes in invasive species distributions and implications for management
All modelled species showed important differences between current and future climate favourability. The variables driving favourability changes measured the severity of both heat-and cold waves (Tmax30 and TmaxAvemin), and the frequency of droughts (Pr1mm5d) and periods of substantial precipitation (PrAvemax5d). However, their contribution was not straightforward to all modelled species, likely because of the differences in their ecology. While both gains and losses are to be expected, losses generally outweigh gains, supporting previous evidence that invasive species will not always benefit from climate change at the regional level (e.g. Bezeng et al., 2017).
For example, N. vison and E. astrild are projected to mostly experience favourability losses, which might impair their performance and increase their vulnerability to other factors, including management (Bellard et al., 2018;Bezeng et al., 2017). Extensive areas of their current distribution will undergo a marked decrease in the frequency of periods of substantial rains, which are favoured by semi-aquatic mammals and sub-tropical birds inhabiting wet habitats (e.g. Ahlers et al., 2015), coupled with more intense heat waves (negative for the American mink according to our models).
Yet, species such as A. amandava, M. monachus, P. krameri and T. scripta are projected to experience more favourable climate conditions across extensive areas of the Iberian Peninsula in the next decades, thus helping further spread and establishment. According to models, they will be favoured more suitable temperature extremesdecreased/increased frequency of cold-and heat waves, respectively (see also Pérez-Santigosa et al., 2008;Shwartz et al., 2009), and more abundant heavy rains, as projected in the east of the peninsula, which are required by aquatic reptiles (T. scripta; Pérez-Santigosa et al., 2008) or bird species native to monsoonal regions (A. amandava; Stiels et al., 2015). Nevertheless, these four species predominantly spread with both intentional and unintentional human assistance but not through natural dispersal, so human-mediated introduction and spread pathways must be controlled to prevent their establishment in areas with potentially more favourable climate conditions in the future. T. scripta is a species of concern at the EU level, while A. amandava, M. monachus and P. krameri are species of national concern, so they are subject to strict restrictions regarding their possession, transport, trade or release into nature, they must be under surveillance, and their populations must be eradicated, or controlled when eradication is not possible .
Invasive species that are less tolerant to high temperatures and drought may be particularly prone to climate change-induced losses within their southern ranges. Models covering the European continent (if extensive data on both species occurrence and extreme climates were available at this scale of extent and resolution) might anticipate an expansion of these species' potential ranges towards the northeast. This is because correlative models, based on empirical species-environment relationships, simply predict changes in the climate space currently occupied by the species. They do not reveal the underlying mechanisms that may either drive species north, or help them cope with climatic changes and persist under suboptimal conditions. Experimental data would be needed for assessing if and how invasive species can adapt to climate change, including climate extremes, but this information is rarely available. Caution is thus required when predicting that invasive species may disappear from unfavourable areas in the future. At best, they may have to adapt to increasingly harsher conditions; at worst, they may capitalize

| Model performance and limitations
The high consensus obtained among predictions from our selected models, with a single PCA axis accounting for 96% to 98% of their variation under both current and future climate scenarios, allows for great confidence in the robustness of provided predictions. Yet, models need to have their predictive ability thoroughly assessed before their forecasts can be used to inform decisions, and this assessment should address different aspects of model performance. Following best-practice standards, we used a structured process to incorporate information from three evaluation metrics that measured performance from different perspectives during model selection Lake et al., 2020). As we used for cross-validation the spatially explicit block method, which splits the modelled area into geographic regions that contain a similar number of occurrence and absence records, overall model performance is a good indicator of how well our developed SDMs would transfer to unsampled regions or perform under future climate conditions in the same region (Lake et al., 2020;Wenger & Olden, 2012). Transferability is especially important for invasive SDMs that aim at projection into incipient or unoccupied habitats under future conditions (Lake et al., 2020;Morey & Venette, 2020;Wenger & Olden, 2012 (Bezeng et al., 2017). Some SDMs based on invasive species distributions at early stages of invasion have proven capable of accurately forecasting the most likely areas to be colonized next (e.g. Barbet-Massin et al., 2018;Pereira et al., 2020;recent validation by Del Moral et al., 2017of Muñoz & Real, 2006 projections for P. krameri spread), although the contrary has also been reported (Briscoe Runquist et al., 2019). In our study, it was not possible to obtain equivalent data on climate extreme variables outside our study region. In any case, to ensure that using only a part of these species' current distributions did not impair our models, we modelled also their global distributions using bioclimatic variables. We showed that our use of regional occurrence data and solely climate extreme variables did not severely truncate the obtained potential distributions: the globally derived patterns were essentially similar within our study region, with visible exceptions regarding only some populations (Appendix S3, Figure S3. Trachemys scripta predicted results (Estrada et al., 2016;Mendes et al., 2020). Even for invasive species, which are introduced artificially and often at separated points in space, future ranges are inevitably linked to current ranges, from where they can only naturally move to areas they can reach from there, and from where they are more likely to move if the local populations are large and locally widespread. Spatial trend surfaces reflect these patterns (Legendre & Legendre, 2012 Spatial autocorrelation (i.e. the fact that nearby observations tend to be more similar than distant ones, rather than independent) stems from both biotic processes such as breeding, interactions and natural dispersal, and from autocorrelation in environmental drivers such as climate and geomorphology (Legendre & Legendre, 2012). Spatial autocorrelation in model residuals is typical and intrinsic of ecological data, as spatially contagious biotic processes drive the presence of spatial structure that is never completely explainable by environmental variables (Gaspard et al., 2019). Nonetheless, it does not invalidate model predictions (Diniz-Filho et al., 2003). As our climatic models did exhibit residual spatial autocorrelation, we used the spatial trend of each species (Appendix S3, Figure S3.2) to incorporate spatial processes not accounted for by the modelled variables (Legendre & Legendre, 2012; see section 2.3.2). Furthermore, we selected the models for the final ensembles using cross-validation on spatial blocks (Appendix S3, Figure S3.1) rather than on random intertwined samples of localities, to avoid the underestimation of prediction error (which is common with random cross-validation in the presence of spatial autocorrelation) and consequently inappropriate model selection (Valavi et al., 2019).
Model predictions might be improved by considering also nonclimatic variables, such as land use, human population density and other anthropogenic factors. However, the aim of the present study was to assess potential distributional changes from the places where these species are currently established and prone to disperse naturally in response to climate extreme events. Land use and other human variables would obscure these effects. Models based solely on climate variables have previously proven able to forecast the spread of invasive species with considerable accuracy (e.g. Barbet-Massin et al., 2018). While the importance of climate-based models should not be underestimated, model predictions should be interpreted with caution and considered as first approximations to the potential magnitude and broad patterns of future impacts, rather than as accurate simulations of future species distributions (Pearson & Dawson, 2003).

| CON CLUS IONS
Understanding how invasive species might change their potential distributional area in the future is critical for the implementation of adequate surveillance, control or eradication measures into local and regional management plans. Fitted SDMs indicate that the selected invasive species can undergo significant changes in their current potential distributions, with losses potentially outweighing gains.
However, four species were projected to experience more favourable climate conditions in extensive areas of the Iberian Peninsula in the future, which should be the focus of intensive surveillance and monitoring. This information is thus useful both to prevent further spread of current invaders, and to manage their populations in areas under decreasing habitat favourability, where they may be easier to target. The methodological framework presented here can also be applied to different taxonomic groups and geographic scopes to inform on areas potentially threaten by biological invasions in a changing climate.

ACK N OWLED G M ENTS
This study was partially funded by the regional Government of within the framework of projects POCI-01-0145-FEDER-030931 and PTDC/BIA-ECO/0207/2020.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The environmental data that support the findings of this study are available in the DRYAD public repository under the following DOI: https://doi.org/10.5061/dryad.rxwdb rv8x. The R code used in our analysis is provided in a public GitHub repository: https://github. com/AMBar bosa/extre me-invaders. The species occurrence data were compiled from multiple different sources (see Appendix S1) and are available from the authors upon request.