Integrating climate and land‐use change scenarios in modelling the future spread of invasive squirrels in Italy

The establishment and spread of invasive alien species may be influenced by several mutually interacting factors, whose understanding is paramount to develop effective biosecurity policies. However, studies focused on modelling spatially explicit patterns of future invasion risk have so far focused on species response to climate change impacts, while land‐use change has been neglected. We investigated how the interplay between climate and land‐use change could affect the future potential distribution and dispersal corridors of four alien squirrels introduced to Europe (Sciurus carolinensis, Callosciurus finlaysonii, Callosciurus erythraeus and Eutamias sibiricus).

Climate and land cover change, as well as their mutual interplay, may promote and increase the establishment and spread of alien species (Dukes & Mooney, 1999). However, studies focused on modelling spatially explicit patterns of future invasion risk have so far mostly relied on climate change alone (Cheung et al., 2009;McDonald & Brown, 1992;Pearson & Dawson, 2003), probably because landuse change scenarios have been less available and reliable than climate change scenarios, at least up to 2000s (e.g., Verburg, Schulp, Witte, & Veldkamp, 2006). While climate change scenarios under different greenhouse emission rates are regularly produced by the Intergovernmental Panel on Climate Change (IPCC; e.g., IPCC, 2007), land-use change scenarios are difficult to produce, considering the range of factors influencing future landscape transformations (but see Rounsevell et al., 2006;Verburg et al., 2006). Therefore, studies taking into account both climate and land-use change are rare and focused on change in biodiversity (Oliver & Morecroft, 2014;Titeux et al., 2017) and only Bellard et al. (2013) used this approach in evaluating suitable areas for introduced species.
This overall lack of integration is particularly concerning because the interactions among multiple drivers of global change are one of the main causes of uncertainty in predicting effects on biodiversity (de Chazal & Rounsevell, 2009;Parmesan et al., 2013). The lack of integrated approaches may hamper a full understanding of how interactions among drivers may affect biodiversity changes, leading to unreliable predictions and misleading conservation recommendations (Sirami et al., 2017;Titeux et al., 2017). Interacting effects of climate and land-use change may accelerate the current dramatic rate of biological invasions, favouring the establishment and spread of introduced species. However, it could also slow down, hamper or revert the invasion by species unable to cope with the new conditions. Among mammal invaders, squirrels and commensal murids represent a major threat worldwide, being highly successful in establishing populations after introduction (Capizzi, Bertolino, & Mortelliti, 2014;Howald et al., 2007;Palmer, Koprowski, & Pernas, 2007). About 250 squirrel introduction events by at least 20 species have occurred in the last 150 years (Bertolino, 2009;, most of which have been successful (Bertolino, 2009). The impacts of alien squirrels range from competitive exclusion of native species (Gurnell, Wauters, Lurz, & Tosi, 2004;Mazzamuto, Morandini, et al., 2017), to tree debarking (Kuo, 1982;Mayle, Proudfoot, & Poole, 2009;Mori, Mazzoglio, Rima, Aloise, & Bertolino, 2016), and parasite/disease transmission (Marsot et al., 2013;Romeo et al., 2014;Tompkins, Sainsbury, Nettleton, Buxton, & Gurnell, 2002).
The squirrel species introduced to Europe are mainly arboreal (Ancillotto, Notomista, Mori, Bertolino, & Russo, 2018;Palmer et al., 2007), and often originate from other bioclimatic realms (Di Febbraro et al., 2013;Di Febbraro, Martinoli, Russo, Preatoni, & Bertolino, 2016). Invaded environments may be differently permeable to native and alien squirrel range expansions, according to the extent and spatial patterns of forests and other key landscape features. Squirrels do not easily move through non-forested matrices (Bakker & Van Vuren, 2004) and their ability in crossing forest gaps or using tree corridors varies across species (Bowman & Fahrig, 2002;Bridgman et al., 2012;Stevenson et al., 2013). Most European populations of introduced squirrels are expanding their distribution range (Bertolino, Cordero di Montezemolo, Preatoni, Wauters, & Martinoli, 2014;Dozières et al., 2015;Goldstein, Butler, & Lawton, 2016;Gurnell, Lurz, & Bertoldi, 2014), offering a suitable model to test the invasion risk they pose considering the effects of both climate and land-use change. Italy makes an ideal setting to test such processes: the region offers important ecological corridors through the vast, continuous stretches of forests that run along the Apennine chain from north to south.
Besides, Italy is home to populations of all alien squirrel species recorded in continental Europe (Bertolino, Colangelo, Mori, & Capizzi, 2015). The Eastern grey squirrel Sciurus carolinensis and the Siberian chipmunk Eutamias sibiricus originate from temperate deciduous forests (respectively, of North America and Asia), with chipmunks being recorded up to the taiga and the southern tundra (Thorington, Koprowski, Steele, & Whatton, 2012 (Thorington et al., 2012), suggesting that they may respond differently to the same environments where the Eastern grey squirrel has also been introduced. The only other alien squirrel species in Europe, the Barbary ground squirrel Atlantoxerus getulus, is present only in Fuerteventura, Canary Islands (Lopez-Darias & Nogales, 2008). All such species are native to the Northern Hemisphere (Thorington et al., 2012). Italy is a long and narrow peninsula with a strong northsouth gradient in environmental parameters, including climate and land cover. This could influence the way the different squirrel species will spread according to their colonization abilities. In a recent paper, Di Febbraro et al. (2016) detected differential responses to climate change of squirrels introduced worldwide, with a pool of species undergoing an increase in their predicted future non-native distribution and others retreating from the invaded ranges.

K E Y W O R D S
Alien species, biological invasions, connectivity models, global change drivers, Sciuridae In this study, we explored how the interplay between climate and land-use change will influence the future potential distribution and dispersal corridors of alien squirrels in Italy. We used Species Distribution Models (hereafter, SDMs; Guisan & Zimmermann, 2000) and circuit theory methods (McRae, Dickson, Keitt, & Shah, 2008) to test whether future scenarios based, respectively, on climate change alone and climate change and land-use scenarios will lead to different effects on the range and connectivity of alien squirrel populations. Specific objectives of this study were the following: (a) to assess the current potential distribution and dispersal corridors for alien squirrels in Italy, considering both climatic and landuse variables, (b) to predict the potential distribution and dispersal corridors of these species under future scenarios based on climate change only and on pooled climate and land-use change, and (c) to quantify differences among species response to both of the considered global change drivers. Accordingly, we predicted (a) a speciesspecific response to interacting climate and land-use change, with some species favoured and others limited by the increase in temperature (IPCC, 2007) and the reduction of forest and semi-natural areas (Verburg et al., 2006), and (b) that effects of interacting climate and land-use change may overturn predictions conducted under climate change only.
Spatial behaviour varies among these species: dispersal distances may vary from 170 m for the Siberian chipmunk (Marmet, Pisanu, & Chapuis, 2011), up to 5 km (mean 1-1.5 km) for the Pallas's squirrel (Guichon & Doncaster, 2008) and up to 7 km for the grey squirrel (Okubo, Maini, Williamson, & Murray, 1989). Grey and Pallas's squirrels may also spread where the suitable forest habitat is highly fragmented, albeit at a lower speed Bridgman et al., 2012). Conversely, the Siberian chipmunk is typical of glades and avoids the agricultural matrix (Jo, Seomun, & Baccus, 2014). All these species are known to adapt also to urban and suburban areas (Bertolino, 2009;Marmet et al., 2011).

| Analytical framework
To assess the effect of interplaying climate and land-use change on the studied species, we modelled their distribution and dispersal corridors under present-day environmental conditions and predicted their alterations under two 2050 climate change scenarios (IPCC, 2007) and two 2030 land-use change scenarios (Verburg et al., 2006).
Since the study area covers just a small portion of the global ranges of the species under analysis, SDMs were calibrated using a hierarchical structure, from global to regional scales (Gallien, Douzet, Pratte, Zimmermann, & Thuiller, 2012;Smeraldo et al., 2017), to avoid truncated niche estimations (Barbet-Massin, Thuiller, & Jiguet, 2010;Raes, 2012). Accordingly, a first set of models was fitted considering the global species range and bioclimatic predictors (i.e., global SDMs, see Supporting Information Appendix S1 for further details). Then, we trained a second set of models, refining predictions at the study area level (Italy, i.e., regional SDMs, Supporting Information Appendix S1).

| Species occurrence data
We calibrated global SDMs with species occurrence data collected from both native and non-native ranges (Broennimann & Guisan, 2008;Di Febbraro et al., 2013;Mainali et al., 2015). Occurrence records were gathered by sampling randomly the IUCN species' range maps (IUCN, 2012) and adding these records to those from the non-native ranges extracted from the "Global Biodiversity Information Facility" database (Strubbe, Jackson, Groombridge, & Matthysen, 2015) and validated by the authors (Supporting Information Appendix S1). As for regional SDMs, we obtained the occurrence records from: (a) datasets compiled by the authors during long-term (>20 years) monitoring projects (see Martinoli et al., 2010;Bertolino et al., 2014), (b) a citizen-science squirrel project ("SaveRedSquirrels": Mori & Menchetti, 2014), and (c) information uploaded by the public to online platforms (www.iNaturalist.org: Table S1). We examined critically the citizenscience data and used only those that were fully reliable (i.e., with photos attached). Although opportunistic records may provide accurate predictions of species distribution (Tiago, Pereira, & Capiñha, 2017), they are often spatially auto-correlated and/or discontinuous (Boitani et al., 2011) due to a generally unknown and unbalanced sampling effort that can vary widely across space (van Strien, Swaay, & Termaat, 2013). Therefore, the initial occurrence dataset including 1,485 records underwent a filtering procedure (see Supporting Information Appendix S2). After filtering, we obtained a dataset of 37 records for C. erythraeus, 89 records for C. finlaysonii, 50 records for S. carolinensis and 11 records for E. sibiricus.

| Environmental variables
We calibrated global SDMs considering, as the initial set of environmental predictors, 19 bioclimatic variables from the Worldclim database (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005; http://www. worldclim.com/current). Bioclimatic predictors were upscaled at a spatial resolution of ca. 50 km, retaining the following six after checking for multicollinearity (VIF ≤ 5; Zuur, Ieno, & Elphick, 2010): mean diurnal range (bio2), isothermality (bio3), mean temperature of wettest quarter (bio8), precipitation seasonality (bio15), precipitation of warmest quarter (bio18) and precipitation of coldest quarter (bio19). As for regional SDMs, we considered the 19 Worldclim predictors, in combination with 18 land-use categories, all rasterized at a resolution of about 1 km. We calculated land-use predictors from the Corine Land Cover 2012 categories, which describe EU land cover in the 2006-2012 time interval at a spatial resolution of 100 m (European Environmental Agency, 2016). Land-use categories were grouped into broader classes (Supporting Information Table S2) according to the ecology of the studied species and following the categories considered in the future land-use change scenarios (Verburg et al., 2006). The four squirrel species live in forests, and show a different ability to move along tree lines, woodlots or shrubby areas (Ancillotto et al., 2018;Guichon & Doncaster, 2008;Mori et al., 2018;Wauters & Gurnell, 1999), thus forests and semi-natural areas are considered areas suitable for dispersal and establishment.
Artificial areas are semi-permeable barriers since the four squirrel species are known to have established populations in urban areas (Ancillotto et al., 2018;Bertolino et al., 2014;Guichon & Doncaster, 2008;Mori et al., 2018). Variables related to agricultural activities (e.g., permanent crops, rice fields, pastures, arable lands) were not considered because generally they represent a matrix not suitable for squirrels. We finally calculated the Euclidean distance from each of the reclassified categories. As similarly done for global SDMs, the final set of predictors was subselected by checking for multicollinearity and included the following eight variables: isothermality (bio3), mean temperature of wettest quarter (bio8), mean temperature of driest quarter (bio9), precipitation of wettest month (bio13), precipitation seasonality (bio15), Euclidean distance from artificial areas, Euclidean distance from semi-natural areas and Euclidean distance from forests.

| Species Distribution Models
For global SDMs, we generated a committee averaging as the main outcome, instead of a "traditional" occurrence probability prediction (Gallien et al., 2012). The committee averaging describes the percentage of agreement on the species presence among various model projections (Thuiller, Lafourcade, Engler, & Araújo, 2009) and was used to weight background points generated to train regional SDMs (Supporting Information Appendix S1).
For regional SDMs calibration, we followed the "ensemble of small models" approach to avoid model overfitting (Breiner, Guisan, Bergamini, & Nobis, 2015;Di Febbraro et al., 2017), a problem arising when few records are available with respect to the number of predictors (Guisan & Zimmermann, 2000). Accordingly, we calibrated for each species a set of models considering all possible combinations of the eight environmental variables taken by two at a time (for a total of 28 combinations per species), and then averaging the results. Regional SDMs were calibrated with an ensemble forecasting approach using the R package biomod2 . We considered the following five modelling algorithms: generalized linear models (GLM), generalized additive models (GAM), generalized boosted models (GBM), random forests (RF) and maximum entropy models (MAXENT). For each species, we randomly placed a set of 10,000 background points over a region identified by all the WWF Terrestrial Ecoregions (Olson et al., 2001) including species records, following the so-called "BAM" framework (Barve et al., 2011). Such framework suggests that, among the factors determining species distribution (i.e., biotic and abiotic), the area that has been accessible to a species through dispersal over a relevant time should guide the definition of the modelling background area (Barve et al., 2011). For invasive species in non-native range, detailed information on the exact introduction locations as well as on the actual long-term dispersal abilities is often lacking or unavailable, seriously hampering an accurate background area delineation. Therefore, we defined the background area in an ecologically relevant way as suggested by Guisan, Petitpierre, Broennimann, Daehler, and Kueffer (2014), that is identifying the set of WWF Ecoregions within which our studied species occur, instead of considering the whole Italy. To evaluate the predictive performance, we randomly split each occurrence dataset into a 70% sample for model calibration and the remaining 30% for model validation, repeating the procedure ten times and averaging the results. Furthermore, we performed a "block" cross-validation (Muscarella et al., 2014): data were split into four geographically non-overlapping bins of equal numbers of occurrences, corresponding to each corner of the entire geographical space. This method has been used to assess model transferability, that is the ability to extrapolate predictions into new areas (Roberts et al., 2017), and to penalize models based on biologically meaningless predictors (Fourcade, Besnard, & Secondi, 2018). The predictive performance of each model was assessed by measuring the area under the receiver operating characteristic curve (AUC: Hanley & McNeil, 1982) and the true skill statistic (TSS: Allouche, Tsoar, & Kadmon, 2006). To avoid using poorly calibrated models, only projections from models with AUC ≥ 0.7 were considered in further analyses (Di Febbraro et al., 2016). Model averaging was performed by weighting the individual model projections by their AUC scores and averaging the result (Marmion, Parviainen, Luoto, Heikkinen, & Thuiller, 2009). We calculated the relative importance of variables from the ensemble model using the functionality provided in the biomod2 package (Jiguet, Barbet-Massin, & Henry, 2010). Lastly, we calculated the spatial autocorrelation in regional SDMs residuals through Moran's I correlograms (Pottier et al., 2013;Di Febbraro et al., 2015; Supporting Information Appendix S3, Figure S3.1). We projected regional SDMs over two climate change scenarios derived by the fourth assessment of IPCC (IPCC, 2007), as well as two land-use change scenarios developed by Verburg et al. (2006). We considered the scenarios by the fourth IPCC assessment instead of the most up-to-dated fifth one (IPCC, 2013), since these lastly released scenarios are available only for climate change, without any analogous for land-use change.
We set regional SDM projections as to rely alternatively either on climate change (OCC) only or on pooled climate and land-use change (CLUC). For OCC scenarios, we considered climate change forecasting at 2050 and kept land-use predictors at current time values (Supporting Information Table S3). For CLUC scenarios, we pooled climate change predictions at 2050 and land-use change predictions at 2030 (Supporting Information Table S3). Since the scenarios produced by Verburg et al. (2006)

| Potential connectivity corridors
Potential connectivity maps for the current time and future scenarios were built using the circuitscape software (McRae et al ., 2008). This program applies principles borrowed from the electrical circuit theory to generate multiple random walk pathways on a set of habitat nodes and a resistance surface, calculating the relative costs of moving through the entire landscape (McRae et al., 2008). The main outcome is a conductance map quantifying the likelihood of a moving subject choosing to cross a cell relative to others available to it (McRae et al., 2008). Using non-validated expert opinion to develop resistance surfaces represents a major weakness of most landscape resistance modelling efforts (Seoane, Bustamante, & Díaz-Delgado, 2005). Thus, we obtained resistance maps by applying a transformation function on the suitability maps by SDMs (Keeley, Beier, Keeley, & Fagan, 2017;Trainor, Walters, Morris, Sexton, & Moody, 2013). Such function defines an inverse relationship between suitability and resistance values allowing different possible shapes (i.e., from linear to negative exponential; Keeley et al., 2017), in which resistance = 0 when suitability = 1 and resistance = 100 when suitability = 0 (Trainor et al., 2013). We tested three alternative suitability transformations according to three function shapes, that is nearly linear, moderately and highly exponential (see Supporting Information Appendix S4 and Figure   S4.1). Subsequently, we binarized the conductance maps to obtain potential connectivity corridors for the current and future times.
In particular, we included all the pixels in the corridor network whose conductance was greater than the mean conductance value in the study area plus one standard deviation (Elliot, Cushman, Macdonald, & Loveridge, 2014).

| Quantification of climate and land-use change effects on species distribution
For each species and scenario, we calculated two metrics of climate and land-use change effect on both range and corridor binary maps: net change (in terms of gain/loss percentage) and geographical stability (expressed as percentage overlap between current and future maps: Franklin et al., 2013). Furthermore, by combining range binary maps, corridor binary maps, and species occurrence records, we calculated the percentage of potential connected range. Specifically, if a binary corridor by circuitscape intersected a group of species records, all the binary presence patches predicted by SDMs and overlapping this corridor were considered as potentially reachable (i.e., connected) from the areas of observed species occurrence (for further details, see Supporting Information Appendix S5 and Figure S5.1). The five metrics were computed for each combination of binarization thresholds and resistance layers (i.e., 12 combinations).
To quantify the different effects of only climate change and pooled climate and land-use change scenarios on each species, we inserted the values of the five metrics from all the 12 combinations of thresholds and resistance layers in a spider plot. For each species/ scenario/combinations, we calculated the area of the polygon in the spider plot, whose vertices represented the values of the five metrics. Then, for each modelled species, we rescaled areas between 0 TA B L E 1 Predictive performance of global and regional SDMs assessed through random splitting of dataset and "block" cross-validation (the area of the polygon obtained by all the metrics at their minimum possible value, i.e., the most detrimental effect) and 1 (the area of the polygon obtained by all the metrics at their maximum possible value, i.e., the most favourable effect). We tested for statistical differences in polygon areas among the four scenarios by a permutational analysis of variance and Tukey HSD post-hoc test.
F I G U R E 1 Suitability maps predicted by regional SDMs under current environmental conditions. Red colours indicate suitable pixels according to the "mean occurrence probability" binarization threshold, while blue colours refer to unsuitable pixels according to this threshold F I G U R E 2 Conductance maps (grey gradient colours) and suitable range (red and dark green colours) for alien squirrels in Italy generated by circuitscape and regional SDMs under current environmental conditions. Red tones refer to suitable patches that overlaps or are spatially linked to the areas of observed species presence (light green buffers) through binary corridors (yellow polygons). Dark green colours indicate the suitable range fraction that is not spatially linked to any area of observed species occurrence. Conductance values and distribution ranges showed in the figure refer to specific resistance layers ("c = 32" suitability transformation, see Supporting Information Appendix S4) and binarization thresholds ("minimum occurrence probability") To conclude, we ranked the five change metrics according to the effect of only climate change and pooled climate and land-use change scenarios among species. In particular, for each species and scenario, we assigned a +1 score to metrics showing an increasing trend (i.e., net change >0%, geographical stability >50% and connected range >0%), and −1 to those showing a negative trend (i.e., net change <0%, geographical stability <50% and connected range <0%). Then, after summing the scores for each metric among the four species, we obtained a final index where +4 indicates an increasing trend for all the species, +2 an increasing trend for three out of four species, and 0 an increasing trend for two out of four species. The same calculation was applied to negative values, indicating decreasing trends.

| Species Distribution Models
Under both random splitting of the dataset and "block" cross-validation, global and regional SDMs reached good to excellent predictive performances (sensu Landis & Koch, 1977;Swets, 1988), showing AUC values >0.8 and TSS values >0.5 for all the species (Table 1).
Suitability maps are shown in Figure 1, whereas potential connectivity values along with the connected range fraction are shown in Figure 2. Climate variables were more important than land-use variables for all the species except E. sibiricus (Supporting Information Figure S6.1). According to MESS results, regional SDMs reported a negligible extrapolation when predicting occurrence probability in future scenarios (see Supporting Information Appendix S7 and Figures S7.1-4).

| Quantification of climate and land-use change effects on species distribution
Only climate change scenarios predicted a ≥50% range increase for all the species except E. sibiricus, which showed a slight increase only under the B1 scenario. Besides, only climate change scenarios predicted a high range stability (>50%) for all species except C. finlaysonii ( Figure 3). Under these scenarios, we forecasted different, yet limited, effects on connectivity. We predicted a slight increase for C. erythraeus and C. finlaysonii, and a decrease for the other two species. Moreover, all species showed a corridor geographical stability >50%, except C. finlaysonii. All species also showed a reduction

| D ISCUSS I ON
This is the first study combining climate and land-use change to predict effects on alien species range and dispersal corridors, confirming the importance of an integrated approach (Sirami et al., 2017). In fact, unlike previous studies (Lamsal, Kumar, Aryal, & Atreya, 2018;Occhipinti-Ambrogi, 2007;Walther et al., 2009) According to our analysis, the Finlayson's squirrel has the lowest invasion potential. This is a tropical species showing the smallest native range among invasive squirrels (i.e., about 12.5 degrees in latitude: Thorington et al., 2012;Boonkhaw, Prayoon, Kanchanasaka, Hayashi, & Tamura, 2017). Its small native niche may account for its limited invasiveness. Conversely, the grey and the Pallas's squirrels occupy a large native range (about 28 and 26 degrees in latitude, respectively: Thorington et al., 2012), explaining the considerably high invasive potential that our integrated scenario revealed. Yet, uncertainty remains about the relationship between invasiveness and native niche/range size for Callosciurus squirrels, as these come from a geographical hot-spot of tree squirrels' diversity (Koprowski & Nandini, 2008), thus their native niche may underrepresent their actual potential due to high interspecific competition (e.g., Guisan et al., 2014). Moreover, the specific status C. erythraeus is also uncertain, as molecular evidence suggests it may rather represent a species complex (Mazzamuto et al., 2016), that is the native range may be significantly smaller. Despite being the most widely distributed species in the native range (over 44 degrees in latitude: Thorington et al., 2012), the Siberian chipmunk seems to have a very limited spread potential in Italy. This may be due to the low dispersal ability of this ground-dwelling species (80-160 m: Marmet et al., 2011) and its high site fidelity (Marmet, Pisanu, & Chapuis, 2009). Most introduced populations of this species are confined to urban parks, where they rely on food provided by humans (Mori et al., 2018). Moreover, all introduced chipmunks derive from captive-bred individuals and not from wild-trapped ones, which may limit their invasive abilities (Mori et al., 2018). The taxonomic status of the Siberian chipmunk is currently under revision, and the taxon is probably a species complex (Obolenskaya et al., 2009;Patterson & Norris, 2016). Genetic and morphometric analyses suggest that the invasive chipmunks in Europe belong to the Korean lineage (E. s. barberi), whose range includes only 6 degrees in latitude (Obolenskaya et al., 2009). The restricted native range would explain the low invasiveness predicted for Italian populations. In light of these outcomes, our prediction (i) has also been verified.
The warmer climate expected for Italy in the future is predicted to promote the spread of species which live at the same latitude of Italy but extending their native range further south (S. carolinensis), or that occur in warmer environments (two Callosciurus species), while a limited range increase is predicted for E. sibiricus which is native of northern latitudes. This effect, however, is opposite for three out of four species when adding land-use changes, which predicts a retraction of forests and thus of potential corridors. The reduced range loss predicted for C. finlaysonii may be related to the large forest cover present in southern Italy, where the main population is established, which will be only limitedly reduced by future land-use changes. This could explain why in our scenarios the distribution of this species is mainly driven by climate change.
Although we cannot rule out the occurrence of a potential niche change and expansion for invasive squirrels, which might lead to invasion patterns differing from those we predicted (Ancillotto, Strubbe, Menchetti, & Mori, 2015;Strubbe et al., 2015), our integrated approach represents an important step forward to define their potential spread and provide a sound base for prioritize management actions. In general, we confirmed that grey squirrel invasion in northern Italy will take longer than that recorded for the British Isles, hampered by the occurrence of unsuitable agricultural matrix in NW Italy (Lurz et al., 2001). Furthermore, the grey squirrel's invasion of the British Isles' has also been facilitated by the spread of the squirrelpox virus, which greatly affects survival of the native red squirrels (McInnes et al, 2006;Rushton, Lurz, Gurnell, & Fuller, 2000), never detected in Italian populations of grey squirrels (Romeo et al., 2018).
Squirrels are commonly traded as pets and they are likely to establish, spread and become invasive outside their natural range, even when populations originate from few founders (Bertolino, F I G U R E 4 Areas of spider plot polygons generated by the five climate and land-use change metrics for the four scenarios. In the box plots, the central line represents the median, while the box comprises the interquartile range (IQR). The upper whisker represents the 3rd quartile +1.5 IQR while the lower whisker indicates the 1st quartile -1.5 IQR. p Values above the dashed lines indicate the statistical significance of the difference in polygon areas between only climate change scenarios (OCC) and pooled climate and land-use change scenarios (CLUC). Conversely, p value above the solid line refers to the statistical significance of the overall difference in polygon areas among all the four scenarios. Statistical significance of differences in polygon areas was assessed through a permutational analysis of variance and Tukey HSD post-hoc test 2009). As a consequence, the recent list of invasive species of Union concern (European Regulation 1143/2014) includes four squirrel species (including the grey squirrel, the Pallas's squirrel and the Siberian chipmunk) that affects native species through competitive exclusion and parasite-mediated competition (e.g., Rushton et al., 2000;Gurnell et al., 2004;Mazzamuto, Bisi, et al., 2017;Mori et al., 2018). Managing invasive squirrel may be challenging, as they are popular and charismatic, for their appealing appearance (Bertolino & Genovesi, 2003); therefore, predicting their potential spread may help researchers and managers to better design control programs. indicates an increasing (decreasing) trend for all the species, +2 (−2) an increasing (decreasing) trend for three out of four species, and 0 an increasing trend for two out of four species in landscape suitability and connectivity could act synergistically with climate changes in speeding up the potential spread of certain species. Therefore, when possible, we suggest that multiple drivers are considered in evaluating the level of potential invasiveness of species under future scenarios. This will increase the accuracy of predictions, with beneficial effects on risk-assessment quality and invasive species management. Future modelling efforts should also include the effect of interspecific interactions, such as predation and competition (Sheehy & Lawton, 2014), at least at the local scale, which might affect the predictions made neglecting such interactions.

ACK N OWLED G EM ENTS
We thank all the citizens who provided us with squirrel occurrences