Effects of climate change and horticultural use on the spread of naturalized alien garden plants in Europe

Climate warming is supposed to enlarge the area climatically suitable to the naturalization of alien garden plants in temperate regions. However, the effects of a changing climate on the spread of naturalized ornamentals have not been evaluated by spatially and temporarily explicit range modelling at larger scales so far. Here, we assess how climate change and the frequency of cultivation interactively determine the spread of 15 ornamental plants over the 21st century in Europe. We coupled species distribution modelling with simulations of demography and dispersal to predict range dynamics of these species in annual steps across a 250 × 250 m raster of the study area. Models were run under four scenarios of climate warming and six levels of cultivation intensity. Cultivation frequency was implemented as size of the area used for planting a species. Although the area climatically suitable to the 15 species increases, on average, the area predicted to be occupied by them in 2090 shrinks under two of the three climate change scenarios. This contradiction obviously arises from dispersal limitations that were pronounced although we assumed that cultivation is spatially adapting to the changing climate. Cultivation frequency had a much stronger effect on species spread than climate change, and this effect was non-linear. The area occupied increased sharply from low to moderate levels of cultivation intensity, but levelled off afterwards. Our simulations suggest that climate warming will not necessarily foster the spread of alien garden plants in Europe over the next decades. However, climatically suitable areas do increase and hence an invasion debt is likely accumulating. Restricting cultivation of species can be effective in preventing species spread, irrespective of how the climate develops. However, for being successful, they depend on high levels of compliance to keep propagule pressure at a low level.


Introduction
The number of species establishing populations outside their native ranges has increased exponentially during the recent century, and this trend is unlikely to slow down in the future (Seebens et al. 2017). Some of these alien species become invasive, i.e. they manage to spread rapidly and maintain large populations over an extended area (Blackburn et al. 2011). These invaders may threaten native biodiversity (Vilà et al. 2011, Blackburn et al. 2014, Bellard et al. 2016 and can have negative effects on ecosystem services and the economy (Vilà et al. 2010, Pimentel 2011. In the case of vascular plants, horticultural trade is the most important pathway of alien species introductions (van Kleunen et al. 2018). Tens of thousands of introduced plant species are cultivated in private and public gardens and green spaces worldwide and represent a huge pool of potential future invaders (Niinemets and Penuelas 2008). In temperate regions, climate change is expected to increase the likelihood that further species from this pool escape cultivation because many garden plants have been introduced from warmer native ranges (Maurel et al. 2016, Dullinger et al. 2017, Haeuser et al. 2018. As human cultivation represents an efficient dispersal pathway, areas becoming suitable under climate change may, moreover, become more rapidly colonized by escaping ornamentals than by non-cultivated species which have to rely on natural dispersal means ( Van der Veken et al. 2008, Corlett andWestcott 2013). However, the implications that cultivation intensity and spatial cultivation patterns may have on the range dynamics of escaping ornamentals under a changing climate have hardly been evaluated so far, at least on a continental extent.
Assessing the relationship between cultivation intensity (i.e. the size and frequency of the populations introduced for horticultural purposes) and spread of alien ornamentals is also topical from the perspective of environmental policy (EU-Regulation 2014). Import or sales bans, voluntary codes of conduct, and raising consumer awareness are the main instruments used to tackle invasions along the horticultural supply chain (Hulme et al. 2017). These measures likely differ in the dampening effect they potentially have on propagule pressure, a factor known to co-determine naturalization success (Lockwood et al. 2005, Feng et al. 2016. What is less well known, however, is whether cultivation intensity is related to the naturalization and spread of ornamental plants in a linear or a non-linear way, which may have considerable implications for the efficiency of 'soft' measures such as codes of conduct or raising consumer awareness. Here, we undertook a large-scale simulation experiment to explore the impact of climate change and cultivation restrictions on the spatio-temporal spread of 15 ornamental plant species (Supplementary material Appendix 1 Table A1) across Europe during the 21st century. We simulated the escape and spatial spread of these species by means of coupled nichedemographic models (Dullinger et al. 2012) run under an assumed constant climate ('baseline') and three scenarios of climate change. The simulation design involved cultivation of the 15 model plants at different levels of intensity which span a gradient from a situation where cultivation restrictions are in force and compliance is high (although never perfect) to an unregulated situation for popular ornamentals. We specifically asked 1) whether species would spread the faster and more successfully the stronger the climate warms, 2) whether increasing cultivation intensity has a linear or a non-linear effect on expansion success and 3) which of the two factors, climate change or cultivation intensity, has the stronger effect on species' expansion.

Material and methods
We simulated the possible 21st century (2010-2090) spread of the selected species across Europe by means of the hybrid model CATS ('cellular automaton-type tool for simulating plant spread') (Dullinger et al. 2012, Hülber et al. 2016. CATS is a cellular-automaton type model that links simulations of demographic and dispersal processes in distinct local populations of a species to the output of species distribution models (SDMs) as indicators of climatic suitability. Spatial spread of the target species starts from areas where garden plants are cultivated and is restricted by both climatic conditions and the pattern of appropriate land cover types. For an overview on the simulation approach see Supplementary material Appendix 1 Fig. A1.

Study region
The study region encompasses the member states of the European Union (except for Cyprus and Malta) plus Switzerland, Liechtenstein, Norway, Albania, Macedonia, Bosnia and Herzegovina, Serbia and Montenegro, and covers a terrestrial surface of 6.2 million km 2 .

Species selection
From the 50 species pre-selected in a collaborative research project on potential future ornamental plant invasions we focused on those 15 species that were 1) successfully cultivated in experimental settings (Liu et al. 2016, Haeuser et al. 2017, Conti et al. 2018) so that we could derive the necessary demographic and dispersal parameters from these experiments, 2) represented by at least 50 occurrence records in the Global Biodiversity Information Facility (GBIF)) and 3) will not lose climatically suitable area in Europe under at least one of three climatic scenarios explored in a previous study (Dullinger et al. 2017; Centaurea americana and Solidago ptarmicoides actually show a minor decrease in suitable area, but were kept to avoid further reduction of the sample of species). The latter criterion was used to concentrate on species relevant for regulations in a warming future. The final list included Amaranthus tricolor, Iris domestica, C. americana, Helianthus debilis, Heliotropium arborescens, Isotoma axiallaris, Lilium formosanum, Lobelia inflata, Monarda punctata, Pennisetum macrourum, Petunia integrifolia, Rudbeckia fulgida, Rudbeckia triloba, S. ptarmicoides and Verbena rigida. All selected species have not naturalized in Europe yet, but are currently cultivated at least somewhere on the continent (Cullen et al. 2011). They have however, naturalized elsewhere in the world, and hence proven their naturalization capacity (Williamson 1999, van Kleunen et al. 2015. The species set includes annual and perennial herbs and grasses which differ in a number of demographic and dispersal related traits as well as in habitat affinities, and originate from different parts of the world (Supplementary material Appendix 1 Table A2).

Species distribution data and climatic maps
We extracted data on the global distributions of the 15 selected species from GBIF (<www.gbif.org/>). Multiple occurrences within 10′ × 10′ grid cells and clearly erroneous records, e.g. in water bodies, were removed. We did not limit records to those from the native range because species are known to partly expand their realized climatic niches in the naturalization range (Petitpierre et al. 2012, Early and Sax 2014, Dullinger et al. 2016. The combination of native and naturalized distributions hence results in the broadest possible representation of the species' realized niche given the global distribution it has reached so far. For characterizing current climatic conditions, we used six bioclimatic variables (data averaged for 1950-2000) provided by WorldClim (Hijmans et al. 2005) at a spatial resolution of 10 geographical minutes: mean diurnal temperature range, minimum temperature of coldest month, mean temperature of warmest quarter, annual precipitation, precipitation of driest month and recipitation seasonality. Pearson's correlations among these variables were <0.7 throughout.
Possible future climates in Europe were represented by three emission scenarios of the IPCC5 report (IPCC 2013) representing mild (RCP 2.6), intermediate (RCP 4.5) and severe (RCP 8.5) climate change. The respective monthly temperature and precipitation time series were taken from the Cordex portal (<http://cordexesg.dmi.dk/esgf-web-fe/ live>; see Supplementary material Appendix 1 Table A3) and used to recalculate 10′ resolution maps of the six bioclimatic variables for future climate scenarios. Fifteen-year running means were then computed to quantify temporal changes in climate over the course of the 21st century at decadal time steps (e.g. 2020 is the 15-year average for 2013-2027, etc.).

Modelling suitable areas
We used the BIOMOD2 framework (Thuiller et al. 2009) to parameterize SDMs by correlating occurrence data from GBIF to the six bioclimatic variables at the 10′ × 10′ resolution. We applied the following modelling algorithms: classification tree analysis (CTA), generalized boosted regression trees (GBM), random forest (RF), generalized linear model (GLM), general additive model (GAM), multiple adaptive regression splines (MARS), flexible discriminant analysis (FDA) and artificial neural network (ANN). Pseudo-absence generation for these modelling algorithms followed bestpractice recommendations by Barbet-Massin et al. (2012) and involved four different approaches: 1) for GLM and GAM we used 10 000 randomly distributed absences; 2) for MARS and FDA 100 randomly distributed absences; 3) for CTA, GBM and RF as many absences as occurrences found in GBIF and selected outside of a radius of 200 km around these occurrences. In the two latter cases, absence generation and hence model calibration, was repeated ten times per species to ensure that selected absences did not bias final predictions; 4) for ANN we used 10 000 absences selected outside of a radius of 200 km around the occurrences and repeated absence generation three times. For all models, the sum of presences was weighted equal to the sum of pseudo-absences. The predictive performance of the models was evaluated by the true skill statistic (TSS; Allouche et al. 2006) based on a three times repeated split-sampling approach in which models were calibrated with 80% of the data and evaluated on the remaining 20%. Evaluated models were then used for ensemble projections (from the techniques combined in each of the four pseudo-absence generation approaches) of the area climatically suitable to each of the 15 plant species under current and future conditions (Araújo and New 2007). Prior to building ensembles, projections of individual models were re-scaled to a common probabilistic scale of [0,1000] and the ensemble projections then computed as weighted averages, with weights determined by the respective TSS scores. Finally, the probabilistic output of these four different ensemble models was aggregated to a total mean. The possibility to combine results of a number of different algorithms into an ensemble projection was also the reason why we preferred binomial SDMs based on pseudo-absences over the use of Poisson point process models (cf. Renner et al. 2015, Guisan et al. 2017. Projections were made for current climatic conditions (year 2010 as baseline climate) as well as for the decadal time steps of the future climates. To produce annual time series of occurrence probabilities from these decadal projections, the 'missing years' in between were linearly interpolated. By using this approach instead of basing SDM projections on annual climatic time series directly we avoided an unduly large effect of annual climatic fluctuations on SDM projections (Hülber et al. 2016).

Species habitat affiliations and European habitat map
We screened online sources (Supplementary material Appendix 1 Table A4) to identify the habitat types known to be suitable to the 15 model species. The suitability of these habitat types was classified using a four-level ordinal scale (with zero meaning unsuitable and three highly suitable) for each species. We cross-tabulated these classifications to the habitat categories distinguished by European CORINE land cover map (CLC, spatial resolution: 100 × 100 m 2 , <www.eea. europa.eu>, re-sampled to 250 × 250 m 2 , cf. Supplementary 4 material Appendix 1 Table A5) and subsequently used CLC to create maps of the distribution of suitable habitat types for each species. These habitat maps were kept constant across the simulation period, i.e. we did not implement any land cover change scenarios into our simulations.

Simulating plant spread
The hybrid model CATS is a tool for simulating shifts in plant species' occurrence and abundance. A detailed description of the model can be found in Dullinger et al. (2012) and Hülber et al. (2016). Briefly, the modelling framework is spatially explicit and discrete in space and time, operating at annual time steps on a two-dimensional raster in which every cell represents an individual site. SDM projections done at the 10′ spatial resolution were resampled to a target grid consisting of 250 × 250 m 2 cells. Each cell of this finer grid was given the suitability value of the 10′ grid it was contained in. Input data for CATS are the initial distributions of the species, climatic suitability (site-specific occurrence probabilities from SDMs), habitat suitability (CLC-based habitat maps), demographic rates and the dispersal matrix. By translating the probabilistic output of the SDMs into demographic rates (such as germination or juvenile survival rates), CATS recomputes the local population structure (i.e. number of seeds in the seedbank, seedlings, juveniles, adults) each year and calculates the populations' annual seed yield. Produced seeds are subsequently dispersed across the raster of sites according to dispersal kernels.

Initial distribution -assumed plant cultivation
We estimated the proportional area generally available for ornamental plant cultivation in each unit of the CLC map on a four-level scale (Dullinger et al. 2017): 0, 0.1, 5 or 10% (Supplementary material Appendix 1 Table A6, Fig. A2). We further assumed that planting of ornamentals in these land cover units is spatially clustered (i.e. not every 250 × 250 m 2 cell of e.g. 'Sport and leisure facilities' has ornamental plants cultivated on 0.1% of its area, but this cultivation is spatially concentrated in 0.1% of the total area of this CLC unit). We implemented this assumption by randomly subsampling from the CLC map equivalent proportions of cells from the respective CLC units for each simulation run. We henceforth refer to this subset of cells as the garden map.
To determine the potential cultivation area of each particular species, we subsequently overlaid the garden map with the SDM projection for the respective species. Because ornamental plants are often cultivated beyond the climatic conditions where they are able to sustain populations in the wild (Van der Veken et al. 2008), climatic suitability was defined in a liberal way when determining the potential cultivation area. Precisely, we used the lowest occurrence probability predicted by ensemble SDMs for any documented presence of the species in the parameterization data (i.e. the GBIF records) as a threshold to delimit the species-specific potential cultivation area on the garden map.
Cultivation frequency levels within the potential cultivation area were derived from the proportion of European nurseries that have particular ornamental plant species on sale. The respective data were taken from Van der Veken et al. (2008) who had analysed 13 000 ornamental plant species and their presence in 250 European nurseries. This database documents a skewed distribution with few species sold in more than 10% of the nurseries, and many available in only one or few nurseries. Based on this empirical distribution we defined six levels of cultivation frequency for our simulation design: species are cultivated in randomly selected cells totalling 0.01, 0.1, 0.5, 1, 5 and 10% of the potential cultivation area. Within these finally selected gardening areas we assumed that individual plants are cultivated at an intensity of 1% of the species' climatically determined local carrying capacity (Supplementary material Appendix 1 Table A7).
Across time, i.e. during the simulation period, gardening habits were allowed to change. In particular, a randomly selected third (for annual plants) or tenth (for perennial plants) of the cells where the species had been cultivated in the previous year was replaced by newly selected sites each year. New cultivation sites were taken only from those parts of the garden map that were climatically suitable to the species according to SDM projections updated for the respective year. Thus, assumed gardening habits followed the changing climate.

Demographic rates
Climatic dependence of local demography of each study species at each 250 × 250 m 2 site was modelled by linking demographic rates (germination, survival, fecundity and clonal reproduction) to occurrence probabilities predicted by SDMs by means of sigmoidal functions (see Supplementary material Appendix 1, and Hülber et al. 2016 for details). Species-specific values of demographic rates were partly derived from results of common garden experiments (Liu et al. 2016, Haeuser et al. 2017) and data gaps were complemented by information from online databases (Supplementary material Appendix 1 Table A2).

Dispersal modelling
We modelled seed dispersal using a compound kernel of four different vectors, namely wind, exo-and endozoochoric dispersal by large herbivores, and human dispersal along streets. Details on dispersal modelling are provided in the Supplementary material Appendix 1, as well as in Dullinger et al. (2015) and Hülber et al. (2016). We partitioned the overall seed yield of the population in each occupied cell among kernels by randomly assigning a fraction between 0.5 and 1% of the seed yield to exozoochoric, endozoochoric and human dispersal, respectively. Remaining seeds were dispersed by the spatially more restricted wind kernel.

Simulation design and statistical analyses of results
For each combination of 15 species, four climate scenarios and six cultivation frequencies we ran ten replicate simulations. For each run, we used the number of 250 × 250 m 2 cells predicted to be occupied by a species in 2090 outside of the garden map as the metric to quantify its invasion success. We compared this metric among climate scenarios and cultivation frequencies by means of linear regression models. Regression models were ran separately with (a) climate scenarios, (b) levels of cultivation frequency, (c) an additive combination of both variables and (d) an additive combination of both variables and their interaction (= the full model) as predictors. For all models we calculated R 2 -values and AICs. Further, to test for non-linearity of the effect of cultivation frequency, we ran the full model with cultivation frequency either logtransformed, square-root transformed or untransformed and compared the alternative models by means of R 2 and AIC values. Prior to the analyses, the number of cells occupied in 2090 was scaled to zero mean and a standard deviation of 1, separately for each species (i.e. for the 240 replicate runs per individual species). This separate scaling filtered speciesspecific determinants of invasion success (demographic and dispersal parameters, habitat affiliations) and allowed to focus on the effects of climate change and cultivation frequency.

Sensitivity analyses
We evaluated the sensitivity of our results to uncertainty in demographic and dispersal parameters as well as in habitat affiliations of species, the garden map and gardening habits by changing the associated parameters within plausible ranges and repeating simulations. Methods and results of these sensitivity analyses are presented in the Supplementary material Appendix 2 Fig. A3-A6 , Table A8. In essence, these analyses demonstrate that the reported results are qualitatively robust against parameter variation.

Species distribution modelling (see Supplementary material
Appendix 3 Table A9 for evaluation statistics) confirms that most species will gain climatically suitable area by the end of the century under at least one scenario of climate change (9, 12 and 8 species under RCP 2.6, RCP 4.5 and RCP 8.5, respectively, see Supplementary material Appendix 3 Fig. A7). However, despite the increase in suitable area, the effect of a changing climate on the spatio-temporal spread of species is ambiguous. Although under the intermediate RCP 4.5 scenario the area occupied during the upcoming decades and, eventually, in 2090 is significantly larger than under the baseline climate, the opposite happens under the mild and severe scenarios RCP 2.6 and 8.5, respectively (Fig. 1A, Table 1). Under all three scenarios, species occupy part of their ranges only transiently (Fig. 1B, Supplementary material Appendix 3 Fig. A9). Loss rates from one decade to next increase sharply at the beginning of the simulation period, obviously a result of fluctuating occupancy of marginally (un)suitable sites and thus pronounced demographic stochasticity. Towards the end of the century, however, loss rates tend to decrease under the RCP 4.5 scenario, while they start increasing again under RCP 2.6 and 8.5 (Fig. 1B). The effect of cultivation frequency on the area occupied by naturalized populations of the 15 species in 2090 is nonlinear, with a log-transformation of cultivation frequency explaining simulation results best (R 2 and AIC values for regression models with linear, square root and log-transformed cultivation frequency were R 2 = 0.45 and AIC = 8028, R 2 = 0.56 and AIC = 7259, and R 2 = 0.62 and AIC = 6737, respectively). Spread rates increase sharply between 0.01% and 1% of cultivation frequency, but this increase levels off subsequently ( Fig. 2A). As compared to climate change, the effect of cultivation frequency on species spread is much stronger. Together, the two factors explain 62% of the variation in the area occupied by wild populations of the 15 species, respectively (Table 1), but the big part of this explained variation is due to the level of cultivation frequency. The interaction between cultivation frequency and climate change scenarios is statistically significant, with restrictions of use having the strongest impact on species expansion under an intermediate climate change scenario (RCP 4.5). However, these interactive effects are much weaker than the main effect of cultivation frequency (Table 1). Across all species, the evaluated three orders of magnitude increase in cultivation frequency resulted in six orders of magnitude increase in the area simulated to be occupied by escaped populations in 2090 (Supplementary material Appendix 3 Fig. A8). This huge effect is mainly driven by three species (Heliotropium  arborescens, Pennisetum macrourum, Verbena rigida), which profit disproportionately from higher cultivation frequency (Supplementary material Appendix 3 Fig. A10). In terms of species richness, the maximum number of species simulated as established in 2090 rises from 7 to 12 per 20 × 20 km cell when increasing cultivation frequency from 0.01 to 10% (Fig. 3).

Effects of climate warming
In summary, our results suggest that the spread of potentially invasive garden plants in Europe will not necessarily be fostered by a changing climate until the end of the 21st century, even if the area climatically suitable to these ornamentals increases. The discrepancy between larger suitable and smaller occupied areas in two of the three scenarios indicates, however, that a warming climate may foster the accumulation of an invasion debt (Essl et al. 2011). The debt arises, first, because ornamental plants do not spread fast enough to track the shifting spatial pattern of climatically suitable sites as suggested by the still strong increase in the number of occupied cells towards the end of the simulation period (Supplementary material Appendix 3 Fig. A8). Second, suitable sites become spatially displaced under a changing climate. This displacement triggers extinction of already naturalized populations in areas becoming unsuitable to the species as indicated by increasing loss rates of once occupied sites towards the end of the century in the mild and severe scenarios RCP 2.6 and 8.5. In the intermediate RCP 4.5 scenario, the effect of these extinctions is apparently outweighed by the larger average increase in the area suitable to the species (Fig. 1A). In the other two scenarios, the net balance is negative, either because the amount of suitable area added is smaller (RCP 2.6) or because the spatial displacement of suitable areas is more pronounced (RCP 8.5).
The accumulation of an invasion debt is in line with the assumption that plant migration will generally lag behind the expected rates of climatic shifts (Loarie et al 2009, Bertrand et al. 2011, Corlett and Westcott 2013, Dullinger et al. 2015. It is more surprising, however, to see these delays arise even though our simulations assumed human gardening habits that adapt to the changing climate and should thus relax dispersal limitations (Van der Veken et al. 2008, Bradley et al. 2012, Corlett and Westcott 2013. While the 'targeted' propagule supply by gardeners certainly facilitates climate tracking, it is hence obviously still insufficient to completely avoid 'climatic debts' (Devictor et al. 2012). These results support findings of other modelling work suggesting that the spread of alien species in a region commonly gains momentum only after a sufficiently large naturalized population has established (Mang et al. 2018). Shifting the spatial pattern of climatically suitable areas implies that this prior build-up of naturalized populations has to start anew repeatedly, especially when the newly suitable regions are separated from already occupied terrain by barriers like mountain ranges, oceans or large tracts of unsuitable habitat. In addition, part of the populations that have already established become less efficient propagule sources for further spread, because the once suitable climate at their sites, or in their surroundings, changes for the worse. This effect may be attenuated by evolutionary adaptation of once-established regional populations to the changing climatic conditions (van Boheemen et al. 2019), but whether and to which degree such rapid evolutionary dynamics will occur is hard to predict.
Our results may underestimate climate change effects on the naturalization and spread of ornamental invasive plants for two reasons. First, the simulations only accounted for changes of the area climatically suitable to the alien species and did not incorporate possible alterations of native vegetation under climate change. In particular, it has been argued that climate warming will cause a pronounced disequilibrium between climatic conditions and vegetation during the next 50-200 years (Svenning and Sandel 2013). An increased frequency of instable, successional stages may increase the invasibility even of those vegetation types that are commonly rather resistant to invaders such as late successional forests (Eschtruth andBattles 2009, Early et al. 2016). Second, the changing climate may be associated with altered land use regimes (Spangenberg et al. 2012) that have also been neglected in our simulations. These land use changes may foster the invasibility of European landscapes, for example where biofuel crop plantations replace former grasslands (Chytrý et al. 2012). However, scenarios of land use change are highly uncertain (Rounsevell et al. 2014), and future development may also increase the resistance of resident vegetation to new invaders, e.g. where economically marginal areas are abandoned (Chytrý et al. 2012).
We emphasize, on the other hand, that our selection of model species excluded those that will lose climatically suitable area in Europe under climate change. While those species are likely fewer in number than the 'winners' (Dullinger et al. 2017), our results are hence not representative for the entire pool of non-native garden plants on the continent. Taking species at random from this pool would have probably resulted in an even stronger shrinkage of the area simulated to be occupied at the end of the century under a changing versus a constant climate.

Effects of cultivation frequency
While our simulations do not support the concern that climate change will markedly accelerate invasion of alien ornamentals in Europe until the end of the century, they demonstrate that restrictions of use have strong effects on the spread of introduced garden plants and that this effect is non-linear. In particular, they suggest that even a relatively low level of cultivation, like the use in 0.1% of the (climatically suitable) gardens, can already trigger the spread of some species quite effectively. Which species become successful invaders under which levels of propagule pressure will often depend on intricate combinations of factors such as climatic suitability, the size and fragmentation of appropriate habitats, the possession of particular traits (Conti et al. 2018), trait complementarity to native species (Carboni et al. 2016), absence of herbivores or pathogens etc. In general, however, our simulation results let it appear unlikely that 'soft' measures such as voluntary codes of conduct alone are able to prevent future invasions of ornamental plants. Rather, they suggest that different policy instruments like risk assessments, legal regulations, industry codes of conduct and raising customer awareness will have to be combined (Hulme et al. 2017) to reduce invasion risk from garden plants markedly. On the other hand, our simulations also emphasize that the effectiveness of such measures is largely independent of how the future climate will develop.

Caveats
As with most models of species range shifts (Zurell et al. 2016), the results of our simulations are subject to parameter uncertainty. First, the SDMs we used may have over-predicted the ranges climatically suitable to the species for two reasons. On the one hand, occurrence points in GBIF did not contain information on whether the population recorded was established/reproducing. We have hence likely included casual populations beyond the climatic limits of self-sustaining populations or even cultivated individuals into model parameterization (Dullinger et al. 2009). On the other hand, fitting models with occurrence data only required assuming absences at sites where no information about the species is 9 available (= pseudo-absences). A number of different strategies for selecting such sites have been developed. Here, we followed best-practice recommendations of Barbet-Massin et al. (2012) in adapting strategies to particular SDM algorithms. However, part of these strategies may lead to over-predictions of suitable ranges (Hanberry et al. 2012). As a corollary, our simulations may also over-rather than underestimate the extent of species' range expansion, independent of how the climate develops. Second, some demographic and dispersal parameters as well as parameters representing habitat affiliation of species and gardening habits were based on rough estimates. The sensitivity analyses show that our qualitative conclusions are hardly affected by varying these parameters within sensible ranges. Nevertheless, the simulated dynamics of individual species do depend on these parameters and are hence associated with considerable uncertainty. Third, our simulations are restricted to only 15 species. The set of alien garden plants already cultivated in Europe and potentially naturalizing in the future is at least two orders of magnitude larger (Dullinger et al. 2017, Haeuser et al. 2018). We do not know how representative the simulation results achieved with our small sample are for this set of species. However, the relative roles of climate warming and cultivation frequency for future ornamental plant spread differ strongly in our simulations, and this difference was robust against subsampling from our 15 species as well as against parameter variation (Supplementary material Appendix 2). We thus conclude that while the magnitude of cultivation and climatic effects certainly varies hugely among species the relative importance of these two factors for future ornamental plant spread in Europe is likely similar across a broader range of species.

Conclusions
Taken together, our simulations suggest that, on average across the 15 garden plants modelled, climate warming is unlikely to foster 21st century spread strongly on a continental scale. However, the weak response to climate warming during the upcoming decades obviously results from a delay of spread. As a consequence of this delay, an invasion debt is accumulating that will likely trigger an important increase in invasion levels once the climate has stabilized again. Restrictions of use may be highly effective in preventing invasive species spread, irrespective of how the climate develops. However, for being successful, they depend on high levels of compliance to keep propagule pressure at low levels.