The future of invasive terrestrial vertebrates in Europe under climate and land-use change

Predicting suitable locations for invasive alien terrestrial vertebrates (IATV) under different scenarios of global change is essential for local and transboundary management aimed to prevent the spread of invasions. Using a spatial modelling approach adapted to invasive species, we identify range-shifts in suitable areas for 15 of the most harmful IATV in Europe, considering future climate and land-use changes. We predict range contractions for seven of these IATV, expansion for four, and inconclusive outputs for the rest. For most Europe, future aggregated distributions show stable or decreasing trends in total IATV richness. Still, specific regions will increase their suitability for additional IATVs, including some protected and last-of-the-wild areas. Our results are informative for early decision-making and long-term strategies to prevent negative effects of IATV. Our approach is based on publicly available data, so predictions can be revised as new data becomes available.


Introduction
Climate and land-use change, biological invasions, and pollution are the main causes of biodiversity loss and species' distribution change (Tylianakis et al 2008, Pyšek et al 2020. Individually, these processes of change may cause a plethora of pervasive impacts that might be exacerbated by their potential synergistic effects (Brook 2008). Particularly, among other potential impacts, climate and land-use change might create disadvantageous environmental conditions (including habitat loss) for native species, but enable the expansion, establishment, and persistence of invasive alien species (IAS) (Fordham et al 2012). Newly advantageous conditions could influence IAS distribution and abundance and broaden their deleterious impacts on native biodiversity, ecosystem functionality, and human economic and social well-being (Stoett et al 2019).
Among IAS, invasive alien terrestrial vertebrates (IATV) pose a particularly significant threat to native species, through predation, competition, hybridization, disease transmission, and habitat alteration, as well as being responsible for damages to human infrastructure, food crops, and forestry (DAISIE 2009).
Prevention and early detection are crucial to avoid the invasion and spread of IATV to new areas, and are also the most cost-effective way to inform rapid responses (Genovesi and Shine 2004). In this vein, detecting areas where ITAV species currently concentrate, or where they may occur but are not detected yet, are valuable outcomes to focus management actions and to support prevention or mitigation (see Polaina et al 2020). On the other hand, in the light of predicted climate and land-use change, foreseeing how these agents of change could modify the distribution of IATV is fundamental to implement longterm strategies for the management of these species (Bellard et al 2016, Gallardo et al 2017, Beaury et al 2020. In the European Union, these directions are included in the regulation on invasive species, which highlights the need to coordinate actions against IAS between member states (European Union 2014). Spatially explicit models, together with the availability of biotic and abiotic data at large scales, allow building a broad forecast of what the future distribution of IATV may look like on a continental European scale.
In this study, we aim to predict changes in suitable areas for IATV according to two of the main agents of Table 1. List of the 15 invasive alien terrestrial vertebrates (IATV) included in the DAISIE's list of '100 of the worst alien species in Europe' (DAISIE 2009). Native ranges information retrieved from CABI Invasive Species Compendium (2020). global change, climate and land-use. Land-use change has been rarely tested in this type of studies, despite its proved relevance in shaping species distributions (Titeux et al 2016). We estimated shifts in environmentally suitable areas (as a surrogate for range shifts) for each of the 15 IATV listed in the '100 of the most invasive alien species in Europe' , one of the available ranking lists comprehensively assessing the impacts of IAS (DAISIE 2009), under different scenarios of climate and land-use change. The focal IATV comprised nine mammals (sika deer, coypu, American mink, raccoon dog, muskrat, raccoon, brown rat, grey squirrel, and Siberian chipmunk), four birds (Canada goose, ruddy duck, rose-ringed parakeet, and African sacred ibis), one amphibian (American bullfrog), and one reptile (pond slider) (table 1).
We used species distribution models (SDMs), which establish relationships between species reported presences and abiotic and biotic factors that are associated with species' survival and establishment in a given area (Soberón 2010). They are a widely used method to produce estimates of species' range shifts, and to identify the current and future hotspots of species invasion (Bellard et al 2013, Gallardo et al 2017, Beaury et al 2020. We adapted SDMs to the specificities of IAS (iSDMs, Gallien et al 2012, Polaina et al 2020 concerning non-equilibrium with the receptor environment and an environmental niche of combined native and invasive conditions. Thus, we built a two-steps hierarchical iSDM that included information from both native and invasive areas, and added weights to the pseudo-absences to correct for non-equilibrium. From iSDM outputs, we calculated present and future IATV richness, and compared their distributions. We also investigated their association with different European habitats, protected and relatively intact areas. Our approach and results provide valuable information on the expansion trends of harmful IATV species, and may assist longterm local and transboundary decision-making in Europe and other regions of the world.

Datasets
We retrieved IATV location records in Europe (table S1.1 (available online at stacks. iop.org/ERL/16/044004/mmedia)) from the Global Biodiversity Information Facility (GBIF.org 2019) to obtain one species presence per 0.25 × 0.25 grid-cell (see filtering process in appendix S1). We selected eight future scenarios of climate and land-use change resulting from four concentration pathways (RCP 2.6, 4.5, 6.0 and 8.5, representing different socioeconomic transformations through time and their associated greenhouse gases emissions) and three time spans (Present, Future 2050 , and Future 2070 ). For bioclimatic predictors, we additionally considered five global circulation models (GCMs) to include uncertainty and to capture model variability (Sanderson et al 2015), which resulted in 40 climatic trajectories (appendix S1).
Land-use variables were retrieved from the Land-Use Harmonization project, which incorporates the interrelation between land-use and climate changes (Hurtt et al 2011; table S1.2). After discarding collinear predictors (VIF ⩾ 0.4; James et al 2013), predictors retained for modelling were fractions of forested and non-forested primary land, non-forested secondary land, urban area, managed pasture, rangeland and cropland. Bioclimatic variables were extracted from the CHELSA database (Karger et al 2017), including annual trends, seasonality, and extremes of temperature and precipitation (three non-collinear predictors; table S1.2). We included other variables of importance for IATV, such as water availability (EEA 2018, Natural Earth 2018), distance to the coast, topography (LP DAAC 2004), and accessibility to major cities (distance-based, Nelson 2008; table S1.2. All datasets and steps followed to obtain derivate variables are available at Polaina et al (2021).

Model fitting
Models included all mentioned non-collinear predictors (VIF ⩽ 4; table S1.2) of present conditions, species presences in Europe, and 5000 random pseudo-absences distributed over the study area (Barbet-Massin et al 2012). We weighted pseudoabsences based on suitability determined by global climatic-only SDMs that included occurrences within the native and invasive range worldwide (Gallien et al 2012; appendix S2).
We fitted SDMs using the BIOMOD2 package in R-software (Thuiller et al 2016, R Core Team 2020) with five algorithms, three alternative sets of random weighted pseudo-absences, and four crossvalidations with 70% of the data for training and 30% for evaluation, which resulted in 60 model outputs per species. Models with high predictive performance (true skill statistical TSS ⩾ 0.7, or 10% of models with the highest TSS value) were combined into a single ensemble model per species, following the committee averaging procedure (Araújo and New 2007). To assess the importance of land-use variables in model fitting and future forecasts, we fitted alternative models excluding these variables.

Model projections
We converted the continuous predictions (between 0 and 1) of each ensemble model based on present conditions into binary maps (0, unsuitable; 1, suitable) using the threshold that maximized the TSS (Araújo and New 2007). To predict Future 2050 and Future 2070 areas of suitability, we used the corresponding climate and land-use variables for those periods. For each period and RCP (eight scenarios), we obtained a single continuous prediction after averaging the five outputs of different climatic GCMs. For binary projections, a grid-cell was suitable if at least three of the five GCMs per RCP scenario predicted a gridcell as being suitable. We compared binary projections of present conditions with the Future 2050 and Future 2070 to estimate the potential expansion/contraction of suitable areas for each species according to the eight scenarios above (Thuiller et al 2016). We approximated the direction of expansion/contraction by measuring the spatial-shift between centroids of the elliptical standard deviations of the binary projections (Tveite 2018). To evaluate the incidence of extrapolating models to non-analogous environmental regions, we calculated multivariate environmental similarity surfaces (MESS) for each future scenario in relation to the environmental conditions where each species is currently present. Then, we assessed the spatial coincidence between environmentally dissimilar areas and predicted suitable areas for IATV (Elith et al 2010).
We defined areas of IATV spread as those where predicted future richness was higher than under present conditions. Predicted richness was calculated as the addition of all individual species' binary maps for a given time span and scenario. We overlaid areas of IATV spread with the main natural habitats in Europe (The Nature Conservancy 2009) to identify those most likely to be affected. We also overlaid our predicted-richness results with protected areas to evaluate the level of already implemented protection that would be available for the additional threat derived from the spread of IATV (UNEP-WCMC and

Model fitting
Outputs from iSDMs showed medium to very good fit (TSS from 0.63 to 0.93) and low uncertainty (CV ensemble models < 0.5 for 13 species; table S3.1). Temperature and precipitation seasonality were the most influential predictors for 14 and 11 species, respectively, followed by other bioclimatic variables. Accessibility was important for six species and the proportion of urban areas for another three (

Model projections
The different scenarios of climate and land-use change yielded a contraction of suitable areas for seven species across periods and RCPs. Four species showed expanded suitable areas, and four had discordant results according to the four RCP's modelled (figures 1, 2, S3.2 and table S3.3). Overall, the greatest changes in size of suitable areas occurred in the 2070 and RCP 8.5 scenarios (the most extreme one). The rose-ringed parakeet and the American bullfrog had the greatest increase in their potential suitable area (12%-68% and 35%-57%, respectively). Conversely, the largest range contractions were predicted for the raccoon dog (27%-50%), the sika deer (26%-48%) and the grey squirrel (30%-48%) (figures 2, S3.2 and table S3.3). The predicted expansions for the brown rat, the American bullfrog and the pond slider were unanimously north-eastwards, whereas the rose-ringed parakeet showed a slight south-eastern drift. Contracting suitable areas had mixed directions (figure 1). Predictions for different future periods (Future 2050 and Future 2070 ) and RCPs (2.6, 4.5, 6.0 and 8.5) mostly agreed in the trends of expansion/contraction and shift direction (figure 1).
Models excluding land-use variables predicted larger suitable areas for ten species under the present scenario and for seven species in the future. General trends of expansion/contraction and shift direction held for most species, except for the coypu, the muskrat and the ruddy duck (figure S4.1). Extrapolations of our main model to non-analogous areas were uncommon within the suitable areas predicted for most species. Still, more than 20% of the Future 2050 predicted suitable areas corresponded to non-analogous conditions for the raccoon, the Siberian chipmunk, the African sacred ibis and the American bullfrog (max non-analogous area = 63%; table S3.4); and for the raccoon dog, the raccoon, the grey squirrel, the Siberian chipmunk, the roseringed parakeet, the African sacred ibis and the American bullfrog in Future 2070 predictions (max. nonanalogous area = 66%; table S3.5). Dissimilarity values measured as MESS were relatively low within suitable areas compared with the rest of Europe (mean MESS within Europe = −1.16 × 10 8 ; mean MESS within suitable areas = −9 × 10 4 ; tables S3.4 and S3.5).
Comparing predicted richness in the present with that in the future, we found that it remained the same for 47%-56% of Europe surface, decreased for 27%-37%, and increased for 13%-23%, according to different climate and land-use change scenarios ( figure S3.3). For example, areas of Belgium and northern Fennoscandia had decreasing IATV richness for all periods and RCPs whereas parts of Italy and Germany had predictions of increased IATV richness for Future 2050 RCP 2.6 and RCP 6.0. The maximum losses and gains of IATV richness corresponded with Future 2070 and RCP 8.5 scenarios. Areas of predicted increase in IATV richness were smaller when land-use was excluded in all scenarios (average decrease of 37%; table S4.4). Areas of potential decline in IATV richness were on average 11% larger than that defined by models including land-use predictors, and areas without changes in richness were ∼3% larger (table S4.4).
Considering predictions per period, summarized as grid-cells where at least two RCPs predicted an increase in the number of IATV, areas of IATV spread mostly overlapped with temperate broadleaf and mixed forest (64% in Future 2050 , 63% in Future 2070 ), followed far behind by Mediterranean (∼15% overlap) and boreal (∼11% overlap) forests (figure 3). Around 28% of areas of IATV spread spatially overlapped with protected areas of Europe, and about 5% overlaid with last-of-the-wild areas (figure 3).

Discussion
It is critical for early decision-making and management prioritization that future IATV expansion or contraction at the individual and aggregated (species richness) levels is understood and accurately predicted. Our focal species provide a good example on how, despite general trends showing general contraction on suitable areas for individual species, the overall IATV richness distribution may still reveal some level of expansion towards new areas. These areas should be particularly monitored and actively prevented from new introductions, especially those placed within protected and relatively intact (last-ofthe-wild) areas.
At the individual species level, previous research showed divergent results for the effects of climate and land-use change on range-shifts, both in vertebrates and plants (Bellard et al 2013, Bezeng et al 2017. Similarly, only four of our studied IATV showed expansion trends in our predictions, with the largest predicted expansions for the rose-ringed parakeet and the American bullfrog. The American bullfrog is one of the main vectors in Europe of the aggressive fungal disease causing the current global decline of amphibians, and the rose-ringed parakeet is a well-known aggressive competitor of native bird species (Menchetti et al 2016, Miaud et al 2016. Consequently, the harmful capacities of each IATV species must be considered individually even in the light of a priori optimistic predictions indicating only modest individual expansions, or even contractions. When applying SDMs for IAS, it is recommended to use global data on native and invasive ranges to approximate the full ecological niche of an invader (Beaumont et al 2009). Yet, the association of IAS with specific land uses and other local factors can differ between native and invasive ranges.  predictors to improve model accuracy, and they were important descriptors of environmental suitability for five IATV.
Identifying changes in IATV richness distribution yields valuable information for decision-making regarding management of biological invasions. Comparing current and future concentrations of IATV species at large scales informs long-term management priorities (Tilman et al 2017). Identified areas of potential increase of IATV richness (figure 3) can focus monitoring campaigns for rapid detection and response. Protected areas where these increases are predicted would be safe if active management actions against invasive species (e.g. restricting landuse changes, intensive monitoring for early IATV detection) are implemented, reducing by around 28% the area where detrimental effects derived from IATV spread are expected. The scarce last-of-the-wild-areas present in Europe would not be notably affected by IATV spread. Areas of France, Belgium, and the Netherlands have a high concentration of IATV species at present, and our predictions show a similar scenario in the future. Therefore, directed management actions would still be necessary to reduce the current and future negative impacts of IATV on native biodiversity.
Errors from model extrapolations to dissimilar environments are moderate for our species. Moreover, invasive species are able to expand into new environmental conditions so dissimilar environments might not be a barrier (Broennimann et al 2007). We acknowledge some unaccounted factors in our research, such as propagule pressure, biotic interactions (e.g. the presence of competitors or facilitators), microclimatic conditions, or differences in evolutionary history, all of which can modify IATV settlement and expansion (Hellmann et al 2008, Yessoufou et al 2014, Ndenda and Yessoufou 2020. Moreover, gaps in species occurrence data still need to be filled, and we need a better understanding of the local features of species ecology within invaded ecosystems (Gethöffer and Siebert 2020, Polaina et al 2020).
Our results provide guidelines for priority setting in transboundary cooperation against relevant IATV. Detecting which species and where they should be firstly tackled at the European level offers an opportunity for managers to act now and implement durable and resilient measures under different climate and land-use scenarios. These guidelines can be continuously updated after incorporating new occurrence data, finer-scale environmental predictors, and expanded to additional species (Beaury et al 2020). We advocate for open-access data to accelerate the incorporation of new information and to assist managers with timely predictions as early actions are vital to combat IAS by using readily available information and tested methods.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.5061/dryad.kkwh70s37.