INTRODUCTION
Studies on plant invasion in temperate mountain regions have traditionally indicated a relatively lower invasion risk in these regions as climatic niche limits are approached with increasing altitude (Alexander et al. Reference Alexander, Naylor, Poll, Edwards and Dietz2009), while human disturbances and propagule pressure decrease (Becker et al. Reference Becker, Dietz, Billeter, Buschmann and Edwards2005; Foxcroft et al. Reference Foxcroft, Jarošík, Pyšek, Richardson and Rouget2010; Haider et al. Reference Haider, Alexander, Dietz, Trepl, Edwards and Kueffer2010, Haider & Kueffer Reference Haider and Kueffer2011). High densities and long residence times of invasive alien plants in the lowlands are correlated with the spread of these species to higher altitudes (Alexander et al. Reference Alexander, Naylor, Poll, Edwards and Dietz2009; Pyšek et al. Reference Pyšek, Jarošík, Pergl and Wild2011) and into protected areas, particularly those exposed to anthropogenic pressures (Pyšek et al. Reference Pyšek, Jarošík and Kučera2002). Several mechanisms may be responsible for this spread, including long-distance dispersal (Pyšek & Prach Reference Pyšek, Prach, de Waal, Child and Brock1994), evolutionary adaptations (Becker et al. Reference Becker, Dietz, Billeter, Buschmann and Edwards2005; Alexander et al. Reference Alexander, Naylor, Poll, Edwards and Dietz2009), low resistance of native communities (Levine et al. Reference Levine, Adler and Yelenik2004), and/or shifts in propagule pressure (Colautti et al. Reference Colautti, Grigorovich and MacIsaac2006), anthropogenic disturbances and climatic regimes (Pauchard et al. Reference Pauchard, Kueffer, Dietz, Daehler, Alexander, Edwards, Arévalo, Cavieres, Guisan, Haider, Jakobs, McDougall, Millar, Naylor, Parks, Rew and Seipel2009; Pyšek et al. Reference Pyšek, Jarošík, Hulme, Kühn, Wilda, Arianoutsou, Bacher, Chiron, Didžiulis, Essl, Genovesi, Gherardi, Hejda, Kark, Lambdon, Desprez-Loustau, Nentwig, Pergl, Poboljšaj, Rabitsch, Roques, Roy, Shirley, Solarz, Vilà and Winter2010).
The Ukrainian Carpathian Mountains (hereafter, UA Carpathians) are a prime example of a temperate mountain range that is experiencing increased levels of invasion by alien plant species established for over a century in the adjacent lowlands and spreading to the interior of the mountains, including existing and proposed conservation areas (Protopopova & Shevera Reference Protopopova and Shevera1998; Protopopova et al. Reference Protopopova, Shevera and Mosyakin2006). The UA Carpathians are rich in biodiversity (Tasienkevych Reference Tasienkevych2008; Prots & Kagalo Reference Prots and Kagalo2012), and over 200 000 ha of protected areas (PAs) have been established (Fig. 1a; Keeton & Crow Reference Keeton, Crow, Soloviy and Keeton2009; Prots et al. Reference Prots, Ivanenko, Yamelynets and Stanciu2010). The Carpathian Convention, established in 2003 to promote sustainable development in the Carpathian Ecoregion, gathered data on biodiversity, landscape features and human development to develop ecological networks in the Romanian, Serbian and UA Carpathians, with the aim to maintain and restore migration corridors between areas of high biodiversity value (Zingstra et al. Reference Zingstra, Seffer, Lasak, Guttova, Baltzer, Bouwma, Walters, Smith, Kitnaes, Predoiu, Prots and Sekulic2009). In the UA Carpathians, the network has been designed to contain the minimum area of overall habitat required to protect predefined high-value conservation features (such as species or interconnected habitats, including PAs, and ecosystem processes) effectively, while being situated in regions that removed from centres of high anthropogenic pressure (Zingstra et al. Reference Zingstra, Seffer, Lasak, Guttova, Baltzer, Bouwma, Walters, Smith, Kitnaes, Predoiu, Prots and Sekulic2009). Over 40% of the UA Carpathians are included in the proposed network, in addition to established PAs (Fig. 1a; Zingstra et al. Reference Zingstra, Seffer, Lasak, Guttova, Baltzer, Bouwma, Walters, Smith, Kitnaes, Predoiu, Prots and Sekulic2009; Deodatus & Protsenko Reference Deodatus and Protsenko2010).
Figure 1 (a) Protected areas/proposed ecological network, showing locations of all 11 invasive study species used for modelling, (b) altitude/hydrology, and (c) roads and settlements in the Ukrainian Carpathians.
Meanwhile, the spread of the invasives outside and inside PAs has been facilitated by anthropogenic disturbances to natural ecosystems, changes to land use in the post-Soviet era and agricultural mismanagement (Török et al. Reference Török, Botta-Dukat, Danczka, Németh, Kiss, Mihály and Magyar2003; Keeton & Crow Reference Keeton, Crow, Soloviy and Keeton2009; Prots & Drescher Reference Prots and Drescher2010; Baumann et al. Reference Baumann, Kuemmerle, Elbakidze, Ozdogan, Radeloff, Keuler, Prishchepov, Kruhlov and Hostert2011). Invasive plant species are found at high densities in close proximity to urban centres, major highways and riparian habitats in the Carpathian foothills and mountain valleys, and disperse into PAs along the frequently disturbed linear corridors, but human settlements and roads, and thus sources of propagule pressure and disturbance (particularly as development and agriculture are not properly regulated by authorities), can also be found within PAs and the ecological network (Fig. 1c). The entire mountain range is well connected through a vast network of roads and rivers (Nazarov et al. Reference Nazarov, Cook and Woodgate2001; see Fig. 1b, c). Plant invasions are likely to intensify in the future, as average monthly temperatures and human infrastructure development in the UA Carpathians as a whole (such as urbanization, road construction and deforestation) and within PAs (foremost tourism) are projected to increase (Nazarov et al. Reference Nazarov, Cook and Woodgate2001; Webster et al. Reference Webster, Holt and Avis2001; Turnock Reference Turnock2002; Bartholy et al. Reference Bartholy, Pongracz, Miklos and Kis2011). An increase in average monthly temperatures may widen the suitable niche space of the invaders (Pauchard et al. Reference Pauchard, Kueffer, Dietz, Daehler, Alexander, Edwards, Arévalo, Cavieres, Guisan, Haider, Jakobs, McDougall, Millar, Naylor, Parks, Rew and Seipel2009), while higher propagule pressure and more frequent anthropogenic disturbance regimes may contribute to the successful introduction and establishment of invasive plants and facilitate their spread, for example by removing natural vegetation and enriching soils (Price Reference Price2006; Pauchard et al. Reference Pauchard, Kueffer, Dietz, Daehler, Alexander, Edwards, Arévalo, Cavieres, Guisan, Haider, Jakobs, McDougall, Millar, Naylor, Parks, Rew and Seipel2009).
However, despite increasing trends of plant invasion, the spatial features, future trends and biodiversity implications of the ongoing invasion process have not been considered sufficiently. Investigating this process is of particular importance because the establishment and subsequent spread of invasive plants in areas of high conservation value has been shown to have negatively affected biodiversity in other mountainous areas (Reinhart et al. Reference Reinhart, Haroldson, Mattson and Gunther2001; Weaver et al. Reference Weaver, Gustafson and Lichthardt2001; Muñoz & Cavieres Reference Muñoz and Cavieres2008), and these impacts may be amplified in the wake of climate change and increases in anthropogenic pressures (Drake et al. Reference Drake, Gonzalez-Meler and Long1997; Dukes & Mooney Reference Dukes and Mooney1999; Dukes Reference Dukes, Mooney and Hobbs2000; Nagel et al. Reference Nagel, Huxman, Griffin and Smith2004; Thuiller et al. Reference Thuiller, Richardson and Midgley2007; Pauchard et al. Reference Pauchard, Kueffer, Dietz, Daehler, Alexander, Edwards, Arévalo, Cavieres, Guisan, Haider, Jakobs, McDougall, Millar, Naylor, Parks, Rew and Seipel2009; Stohlgren et al. Reference Stohlgren, Pyšek, Kartesz, Nishino, Pauchard, Winter, Pino, Richardson, Wilson, Murray, Phillips, Ming-yang, Celesti-Grapow and Font2011).
In this paper, we use Maxent, a species distribution modelling software package designed to work with presence-only data (Elith et al. Reference Elith, Phillips, Hastie, Dudik, Chee and Yates2011), to predict future habitat suitability within the UA Carpathians for 11 widespread and potentially harmful (in terms of loss of local floral diversity; Botta-Dukàt & Balogh Reference Botta-Dukàt and Balogh2008) alien invasive plant species. The model uses six predictor variables, representing gradients in climate, topography and disturbance/propagule pressure to produce projections for future change scenarios, assuming seasonal temperature increases of ≤ 1.8 °C and ≤ 3.8 °C, and increases in anthropogenic disturbances of 10% and 30% for the years 2050 and 2100, respectively.
We hypothesized that, under current environmental conditions, conservation areas were susceptible to invasion by alien plants and that habitat suitability within these areas would increase dramatically in the future. Given that climate, topographic and anthropogenic predictors limit the distribution of suitable habitats for species establishment, we sought to assess the potential intensity (in terms of species and area) of invasion in areas of high conservation importance for both current predictor values and values modified to approximate future climate change and intensification of anthropogenic disturbance.
METHODS
Study area
The study area encompasses the entire range of the UA Carpathians (48°32ʹ N, 23°38ʹ E), which extend over an area of 24 000 km2. The study area lies at an altitude of 95–2030 m, although 94% of the mountains are < 1200 m (Fig. 1b). The highest elevations are located in the southern parts of the UA Carpathians, while the south-west (bordering Romania), west (bordering the Transcarpathian Lowland of Ukraine) and north-west (bordering Poland) Carpathians are characterized by extensive valley systems and relatively gentle slopes. Precipitation of 500–1400 mm yr−1 feeds a dense network of rivers (Fig 1b; Holubets et al. Reference Holubets, Honchar, Komendar, Kutsheraviy and Odynak1988). The July (warmest month) temperature varies from 20° C at the southern edge of the Carpathians and 18° C in the north to 6° C on the highest peaks (Herenchuk Reference Herenchuk1968; Kuemmerle et al. Reference Kuemmerle, Chaskovskyy, Knorn, Radeloff, Kruhlov, Keeton and Hostert2009). Winter temperatures range from –3° C to –10° C. The mountains are dominated by Fagus sylvatica and Picea abies forests, replaced by Pinus mugo and Juniperus communis in the subalpine and grasslands in the alpine belts (Herenchuk Reference Herenchuk1968; Kuemmerle et al. Reference Kuemmerle, Chaskovskyy, Knorn, Radeloff, Kruhlov, Keeton and Hostert2009).
Study species
Using the Alien Plant Ranking System (APRS) developed by the Northern Prairie Wildlife Research Center (APRS Implementation Team 2000) and taking into account reports by regional experts (Botta-Dukàt & Balogh Reference Botta-Dukàt and Balogh2008), we determined the 11 potentially most harmful (to local biodiversity) alien invasive plant species to be Acer negundo L., Ambrosia artemisiifolia L., Echinocystis lobata (Michx.) Torr. & Grey, Helianthus tuberosus L., Heracleum sosnowskyi Manden, Impatiens glandulifera Royle, Reynoutria japonica Houtt., Reynoutria × bohemica Chrtek. & Chrtková, Robinia pseudoacacia L., Solidago canadensis L. and Solidago gigantea Aiton. All species, with the exception of A. artemisiifolia and H. sosnowskyi, were intentionally introduced to Central and Eastern Europe in the 18th and 19th centuries, eventually escaped controlled cultivation, and established self-replicating populations by the beginning of the 20th century (Botta-Dukàt & Balogh Reference Botta-Dukàt and Balogh2008). The species share similar physiological and life history traits (Appendix 1, Table S1, see supplementary material at Journals.cambridge.org/ENC; Walter et al. Reference Walter, Essl, Englisch and Kiehn2005; Weber & Jacobs Reference Weber and Jacobs2005; Botta-Dukàt & Balogh Reference Botta-Dukàt and Balogh2008; Kabuce & Priede Reference Kabuce and Priede2010). All species exhibit fast growth and relatively high reproductive rates. In addition, all species have wide ecological niches and successfully use water and anthropogenic vectors for long-distance dispersal. All the invasives, with the exception of A. artemisiifolia, E. lobata and H. sosnowskyi, reproduce vegetatively.
Combinations of these traits may contribute to competitive dominance of the invasives particularly in early successional habitats (Sakai et al. Reference Sakai, Allendorf, Holt, Lodge, Molofsky, With, Baughman, Cabin, Cohen, Ellstrand, McCauley, O'Neil, Parker, Thompson and Weller2001; Pyšek & Richardson Reference Pyšek, Richardson and Nentwig2007), and populations of the species have been shown to negatively affect ecosystem composition and/or function in temperate invaded habitats outside the UA Carpathians (Botta-Dukàt & Balogh Reference Botta-Dukàt and Balogh2008). For example, A. negundo can spread rapidly in disturbed riparian habitats and is known to prevent regeneration of poplar and willow communities (Mędrzycki Reference Mędrzycki2007). Meanwhile, dense stands of I. glandulifera, another invader of riparian habitats, tend to destabilize riverbanks as the shallow roots of the plants (10–15 cm) do not hold soil efficiently, thus altering abiotic conditions (Hejda & Pyšek Reference Hejda and Pyšek2006; Prots & Drescher Reference Prots and Drescher2010). It must be noted, however, that the impacts of invasive species are largely unquantified in the UA Carpathians and thus require further research.
Data on species’ presence
Several hundred presence records for each of the selected species (Table 1; Fig. 1a) were collected in the field, reviewed in the literature, and supplemented with reliable herbaria records from the University of Lviv (LW), the University of Uzhgorod (UU), the State Museum of Natural History in Lviv (LWS), and the University of Chernivtsi (CHER). Each of the species was modelled individually, with the exception of Solidago spp. and Reynoutria spp. The study species within each of these two genera are very similar in physiology and have overlapping niches in the study region (B. Prots, personal observation; see also Balogh Reference Balogh, Botta-Dukàt and Balogh2008; Botta-Dukàt & Dancza Reference Botta-Dukàt, Dancza, Botta-Dukàt and Balogh2008), and it is thus likely that mistakes exist in older datasets (literature data and field vegetation records) in identifying the genera as separate species. We therefore modelled these genera as one complex.
Table 1 Maxent model results for 11 invasive plant species in the Ukrainian Carpathians - Model performance (AUC score) and relative importance of predictors (in per cent). Potentially highly invasive species, number of presence locations used for Maxent modeling, average model training and test AUCs, and permutation importance for the six predictor variables used for model fitting (in per cent; higher values indicate greater decrease in model performance if predictor is randomly permutated) are shown. Maxtwarm = 40-year average maximum temperature (°C) of the warmest month; mintcold = 40-year average minimum temperature (°C) of the coldest month; s_dist_sett_r - proximity (m) to roads and settlements; s_dist_water = proximity (m) to water bodies; slope = slope (°); sat = yearly sum of daily average active temperatures (°C) above 10° C.
All data were available in a presence-only format. For each species/genus, the majority (> 80 %) of records originated from georeferenced (precision ≥ 10 m) field samples with an approximate mean distance between neighbouring locations of > 1000 m. Of these, 90 % were collected along major environmental gradients relevant to the current distribution of the species (that is, climate and location in relation to water and anthropogenic structures) in order to prevent sampling bias (see Phillips et al. Reference Phillips, Dudik, Elith, Graham, Lehmann, Leathwick and Ferrier2009). The remaining samples consisted of locations identified in herbaria records and confirmed through field observation. Only locations where permanent populations have become established (occupying a sampling unit of 50 m2 in consecutive years) were included in the modelling in order to minimize model inaccuracies due to casual opportunistic observations.
Distribution modelling: Maxent
Maxent modelling is a general-purpose machine-learning method for making inferences from incomplete information. The application has specifically been developed for presence-only data because it does not make assumptions about absences and has been shown to outperform other presence-only modelling methods (Phillips et al. Reference Phillips, Anderson and Schapire2006; Dudik et al. Reference Dudik, Phillips and Schapire2007; Franklin Reference Franklin2009; Elith et al. Reference Elith, Phillips, Hastie, Dudik, Chee and Yates2011). Given a set of grid maps where each pixel represents a local value of a predictor variable and a set of coordinate points depicting species presence, Maxent estimates two probability distributions, one of predictor variables (z) over presence locations, f1(z), and another of predictor variables across randomly chosen background points from the entire study area, f(z), and then determines the values for the response variable (here, suitability of pixels within the study area for the establishment of an invasive plant species) by finding the most uniform distribution of suitable areas given the constraint that the expected value of each predictor under this distribution matches its empirical average at the set of presence locations. The rules, or functions, of how the probability distributions are determined and matched in multivariate space are described by the features, or linear transformations, of potentially complex relationships between the density of predictors and presence/background locations (Elith et al. Reference Elith, Phillips, Hastie, Dudik, Chee and Yates2011).
For each of the nine species/genera, 10 000 random background points were extracted from the study area, representing potential habitat. Due to the geomorphology and presence of extensive river and road systems within the study area, it was assumed that significant geographical barriers to the potential distribution of the invasives do not exist. Along with the background points, 80 % of the presence records were used for model fitting and 20 % for testing. Because the performance of the models is influenced by the particular partitioning step the software assigns to the data, this effect was minimized by running a 5k cross-validation. This method randomly divides the occurrence data into five equal-sized folds, and models are created leaving out each fold in turn. The omitted fold is used for evaluation. A final model run was made for each species using all the presence records for model fitting in order to derive the most robust classification for visual interpretation (Hernandez et al. Reference Hernandez, Graham, Master and Albert2006).
The AUC (area under the curve) statistic obtained by the receiver operating curve (ROC) was used to evaluate the performance of models (Phillips et al. Reference Phillips, Anderson and Schapire2006; see also Evangelista et al. Reference Evangelista, Sunil, Stohlgren, Jarnevich, Crall, Norman and Barnett2008). The ROC plots model sensitivity on the y axis against (1 – specificity) on the x axis for all possible thresholds. Sensitivity is the fraction of presence locations correctly predicted to overlay suitable habitat, and specificity is the fraction of background locations correctly predicted to overlay unsuitable habitat (Fielding & Bell Reference Fielding and Bell1997). An AUC value of 0.9 indicates that 90% of the time when a presence and background location are drawn at random, the first will have a higher predicted suitability value than the second. The statistical significance of the AUC can be determined by comparing the results with random predictions, which would have an AUC of 0.5. Guisan et al. (2007) proposed a classification scheme to assess the significance of AUC values above 0.5, where AUC > 0.90 = excellent, 0.90 > AUC > 0.80 = good, 0.80 > AUC > 0.70 = useful, and AUC < 0.70 = poor (see also Swets Reference Swets1988; Jeschke & Strayer Reference Jeschke and Strayer2008). However, because the AUC does not consider the significance of predicted probability values (Lobo et al. Reference Lobo, Jimenez-Valverde and Real2007), a Wilcoxon ranked sum test implemented with the stats package in the R statistical software (R Development Core Team 2011) was applied for each model to test whether the suitability predictions over presence locations had a higher score than a set of background predictions randomly sampled from the study area (Phillips et al. Reference Phillips, Anderson and Schapire2006).
Maxent also provides a permutation test to assess the relative importance of predictor variables. After a model had been calibrated using model-specific measures of variable contribution (expressed as coefficients), each predictor was in turn randomly permutated at the training points (presences and background) and the decreases in model performance (AUC) were recorded. Permutation values were normalized over all predictors (Phillips Reference Phillips2010; Table 1).
Predictor variables
Initially, 20 bioclimatic and three topographic variables were available for modelling. Of the 20 bioclimatic variables, 19 were retrieved as ESRI grids from the WorldClim global database at a resolution of 1 km (Hijmans et al. Reference Hijmans, Cameron, Parra, Jones and Jarvis2005). The bioclimatic data were regional averages of climatic grids generated by thin spline interpolation of average (1960–1990) monthly climate data from global weather stations. The 20th bioclimatic variable, sum of active temperatures (annual sum of average daily temperatures > 10° C; see Herenchuk Reference Herenchuk1968; Prots & Kagalo Reference Prots and Kagalo2012) was interpolated, using the topo-to-raster function in ArcMap 10, from an ecoregion map of the UA Carpathians (original scale 1:200 000) that displayed climatic regimes within topographic zones. The three topographic variables were derived from vector maps (original scale 1:200 000) on hydrology, roads and settlements (Fig. 1b, c) and from a digital elevation model (DEM; original resolution of 30 m). The vector maps and DEM were provided by the Geography Department of the University of Lviv (Jarvis et al. Reference Jarvis, Reuter, Nelson and Guevara2006; Hostert et al. Reference Hostert, Chaskovskyy, Knorn and Kuemmerle2008; Kruhlov Reference Kruhlov2008; Kuemmerle et al. Reference Kuemmerle, Chaskovskyy, Knorn, Radeloff, Kruhlov, Keeton and Hostert2009; Deodatus & Protsenko Reference Deodatus and Protsenko2010). Layers relevant to propagule pressure and disturbance, such as proximity to water and to settlements and roads, were derived from the hydrology map and the combined map on roads and settlements using the simple (Euclidean) distance function in ArcMap 10. Slope was derived from the DEM.
In order to maintain the spatial accuracy of the DEM, all other variables were rasterized or resampled (in the case of the WorldClim datasets) to a resolution of 30 m using the cubic resampling function in ArcMap 10. The resampling did not improve the resolution of the climatic datasets and was performed solely in order to create raster layers of the same size without generalizing, and thus losing, some of the spatial information provided by the layer with the finest resolution (Yates et al. Reference Yates, McNeill, Elith and Midgley2010). All layers were projected onto the UTM grid, zone 34 with WGS84 datum.
Initial models were run with all variables, and consistently showed that removing irrelevant predictor variables significantly improved the performance of the models. Predictor selection followed three steps.
(1) Initial model: a model was run including all available predictors and the AUC values on training and test data were recorded.
(2) Ecologically-based predictor selection: predictors that approximated (i) limiting factors controlling the ecophysiology of the species and (ii) natural or human-induced disturbances were preferentially selected. For example, because invasive plant species were reported to be limited by extreme temperature regimes (Botta-Dukàt & Balogh Reference Botta-Dukàt and Balogh2008) and because precipitation was not a limiting factor in the UA Carpathians (Herenchuk Reference Herenchuk1968), bioclimatic variables depicting seasonal maxima or minima were selected over annual temperature means and over precipitation in general. Models were rerun with the limited set of predictors, and AUC scores were compared with those derived from the initial model.
(3) Permutation-test based predictor selection: predictors with the lowest permutation importance (< 2 %) in step (2) were removed. A final run was then made and AUC scores were compared with those in previous steps. The scores showed that elimination of irrelevant variables across species had improved the model.
The final predictors were thus: minimum temperature of coldest month (mintcold) as a proxy for susceptibility to frost, maximum temperature of warmest month (maxtwarm) as a proxy for susceptibility to drought, sum of active temperatures > 10°C (sat) as a proxy for length of the growing season, proximity to water (s_dist_water) as a proxy for soil moisture and natural disturbances/propagule pressure, proximity to settlements and roads (s_dist_sett_r) as a proxy for anthropogenic disturbances/propagule pressure, and slope (slope) as a proxy for topographic preference (Table 1; Fig. S1, Appendix 1, see supplementary material at Journals.cambridge.org/ENC, for an example of a Maxent probability model fitted to the six predictors).
Climatic and land-use projections
We applied two simple projection scenarios based on temperature shifts by 2050 and 2100 and high and low anthropogenic pressures. Climate projections assumed the A1B scenario developed by the IPCC Special Report on Emission Scenarios (SRES): an estimated increase in CO2 concentration levels of 532 ppm and 717 ppm by 2050 and 2100, respectively. Based on this scenario, the European ENSEMBLES project developed a series of regional climate models for Hungary (for details see Bartholy et al. Reference Bartholy, Pongracz, Torma, Pieczka, Kardos and Hunyady2009, Reference Bartholy, Pongracz, Torma, Pieczka, Kardos and Hunyady2011). Their calculations extended into the UA Carpathians and were used to adjust the bioclimatic predictor variables according to the proposed 30-year average increases in temperature: 1.8° C and 3.8 ° C in winter and 1.5 ° C and 3.5 ° C in summer for the periods 2021–2050 and 2071–2100, respectively. Based on these increases, new values for the bioclimatic variables were created through simple addition of the mean increases for each pixel value. For example, 18 and 38 were added to all values of mintcold (in ° C × 10) to create the new set of mintcold predictors used for projections for 2050 and 2100, respectively.
In addition, two simple future settings of anthropogenic pressure were developed: (1) disturbances along roads and settlements will not increase above the current level due to low economic development; and (2) more land around settlements and roads will be disturbed due to high economic development (for models of land-use change in Europe, see Rounsevell et al. Reference Rounsevell, Reginster, Araujo, Carter, Dendoncker, Ewert, Hause, Kankaanpa, Leemans, Metzger, Schmit, Smith and Tuck2006; Kuemmerle et al. Reference Kuemmerle, Hostert, Radeloff, van der Linden, Perzanowski and Kruhlov2008). Low and high economic development could also be linked to stronger and weaker nature protection, respectively. The former disturbance setting determined the model projection scenario CL: changes in the climatic regimes for 2050 and 2100 were modelled without incorporating changes to the variables approximating anthropogenic disturbances/propagule pressure. The latter disturbance setting determined the model projection scenario CL&HED: climatic changes were combined with a net decrease in the distance to any given potential human disturbance point because there would be more such points and the impact of existing points may be greater (for a similar approach see Rouget & Richardson Reference Rouget and Richardson2003). That is, for any pixel in the study area, the distance to roads and settlements that had been determined for current conditions decreased by 10 and 30 % by 2050 and 2100, respectively. A paired Wilcoxon signed-rank test (in R) was applied to test whether there was a significant range expansion across species under projections.
The four projections (CL and CL&HED by 2050 and 2100) developed for this study are purely illustrative and are primarily intended to portray general trends in the potential of the study species to profit from climate and land-use changes.
Distribution of suitable habitat in protected areas and the ecological network
In order to quantify the impacts of invasion in PAs and the ecological network in terms of proportion of total area, or pixels, potentially suitable predictions for current distributions and future projections calibrated on all presence points were transformed into binary suitable (= 1) and unsuitable (= 0) values using an optimized threshold based on the ROC curve which maximized sensitivity plus specificity. This approach has been shown to perform well, but is sensitive to low prevalence of occurrence data, which is assumed to be the case here (Liu et al. Reference Liu, Berry, Dawson and Pearson2005). Since information about prevalence could not be derived due to lack of absence data, the results must be considered as relatively conservative binary estimates. Comparisons between species were accomplished by overlaying the binary predictions. Areas that are suitable for one or more species within PAs, thus indicating the range of invasion risk and potential impacts on biodiversity, could then be quantified for current conditions. For future projections, the proposed ecological network was assumed to be a part of areas of high conservation value alongside current PAs.
RESULTS
With AUC values consistently > 0.9 (Table 1) for both training and test data, the predictive performance of the models was excellent. All predictions of habitat suitability were statistically significant (W > 5 000 000, p < 0.001). All predictor variables contributed significantly to model performance, but the relative importance of predictors varied across species and could be explained by species-specific niche preferences.
Binary predictions of habitat suitability clearly suggested a spatial aggregation of suitable habitats for establishment of the species within the study area (Fig. 2). Of the pixels predicted as suitable for the establishment of at least one species in PAs (8 % of the total number of pixels comprising the area within PAs), 13 % overlapped for all nine species (Fig. 2).
Figure 2 Spatial distribution of suitable habitats for potential establishment of at least one or nine invasive plants/genera within the study area and within protected areas under current environmental conditions as determined by the Maxent models.
The spatial distribution of sites suitable for establishment of all nine species/genera showed a preference for major linear habitats along rivers and roads that are in close proximity to centres of high anthropogenic pressure in moist and warm lowlands (up to c. 700 m altitude), namely in the south-west (Upper Tysa Depression and Cirocha-Rika Low Mountains), west (Vyhorlat-Hutyn Volcanic Ridge) and, to a lesser degree, the east (Marginal Beskydy) and south-east (Marginal Gorgany) of the mountains (Fig. 3; Table 2). Suitability decreased at higher altitudes, and PAs at altitudes > 700 m, characterized by relatively steeper slopes and exposed to fewer anthropogenic disturbances, were at lesser risk of being invaded. Furthermore, predictions did not extrapolate far beyond areas in which one or more of the study species were already established (compare Fig. 1 and Fig. 2).
Table 2 Association between the prevalence of invasive plants (average number of all nine study species/genera per pixel) and mean values of maximum temperature of the warmest month (maxtwarm), distance to settlements and roads (s_dist_sett_r), and slope within each of the 33 geomorphologic regions in the Ukrainian Carpathians, listed in descending order by number of all nine study species/genera per pixel.
Figure 3 Biogeographic regions within the Ukrainian Carpathians (see Table 2). Regions shaded grey contain the highest proportion of the nine invasive plant species/genera per pixel.
Under the two scenarios, species may lose climatic and physical/environmental barriers that limit rapid spread today (Table 3). A strong and significant (V = 0, p < 0.001) increase in areas suitable for establishment across species was projected. The net gain of novel suitable habitats in the study area was significantly higher (V = 0, p = 0.002) under the scenario CL&HED than under the scenario CL for both 2050 and 2100. Decreasing the distance to points of anthropogenic pressure by 10% and 30% significantly increased the proportion of the total study area projected as suitable for species establishment (Table 3; Fig. 4). By 2100 (scenario CL&HED), all species were expected to find suitable habitat in at least 20% of conservation areas, and for some invasives significantly greater area became available (for example R. pseudoacacia = 43.3%, H. tuberosus = 30.0%).
Table 3 Current modelled and future projected suitable habitats for invasive species within the Ukrainian Carpathians. Proportions of 30 × 30 m grid cells that contain suitable habitats in the entire study area (SA) and within protected areas (PA) are given for current predictions and future scenarios (CL = climate change; HED = high economic development). All projected models showed significant increases (p < 0.001) in suitable habitats for establishment of invasive species compared to predicted current suitable habitats.
Figure 4 Projected spatial distribution of suitable habitats for establishment of at least one or nine invasive plant species/genera by (a, b) 2050 and (c, d) 2100 within the study area and within protected areas and ecological network assuming (a, c) climate change (CL) and (b, d) climate change and high economic development (CL&HED).
The interior of the UA Carpathians and PAs and the proposed ecological network at higher altitudes became increasingly suitable under a warming climate (Fig. 4). The ecological network may become susceptible to invasion as it extends far beyond current PAs and encompasses areas in the south, west, east and south-east that already function as aggregation centres for a large number of invasives (Figs 3 and 4).
By 2050, when solely climatic changes are expected, 25 % of the total conservation area (PAs and ecological network) was projected to be suitable for at least one species, an increase by 17 percentage points compared with predictions for current conditions. Of this area, 25% became suitable for the establishment of all nine species/genera. Including high economic development and thus increased anthropogenic pressure increased the projected range of establishment of at least one species slightly to 26%, but did not produce greater aggregation of suitable habitats for all nine species/genera (Fig. 4b). The majority of potential habitat.gained was in linear areas in the interior of the mountain range (namely in climatic zones that are currently unsuitable for the establishment of permanent populations). In areas that are highly suitable for invasion today, future projections suggested a lateral spread of species away from major linear habitats and along small waterways (Fig. 4a, b). Projections for 2100 demonstrated a continuation of this trend (Fig. 4c, d).
By 2100, suitable habitats had developed far into the interior of the mountain range, further increasing invasion risk in the central Carpathians and the south (for example in parts of Chyvchyny Crystalline Polonyny; see Table 2). Modelling climate change alone, 44 % of the total conservation area was projected to be suitable for at least one species, and 29 % of these areas were suitable for all nine species/genera. Under the assumption of high economic development by 2100, 47 % of ecologically valuable habitats were projected to be at risk of establishment of at least one species; of these at-risk habitats, 31% were suitable for all nine study species/genera (Fig. 4d). Major watersheds and roads in virtually all PAs and ecological networks were suitable for all nine species/genera; individual species may extend their ranges much further and invade a large number of smaller linear habitats. If current trends continue, by 2100 (under both scenarios), only the most remote parts of PAs at high elevations (such as parts of biogeographic regions Chornohora Polonyny and Internal Gorgany; see Fig. 3) were projected to remain free of invasion.
A comparison between the average habitat suitability gain across species within conservation areas (PAs/ecological network) and the study area as a whole revealed that a relatively greater gain in novel habitats was projected within conservation areas for the future change scenarios (Fig. 5), although, for each species, relatively more pixels were projected to be suitable in the study area compared to conservation areas (Table 3). By 2100, under scenario CL&HED, conservation areas gained 24 % in potentially suitable habitats, as opposed to 22 % in the whole study area. For the same scenario meanwhile, 27 % of the area within conservation areas was on average projected to be suitable, five percentage points less than within the entire study area. A close examination of the spatial patterns of projected potential invasion revealed that the design of ecological networks at low to medium elevations (namely in regions of extensive projected gain in suitable habitats) was responsible for the relatively greater increases in potential invasion risk within conservation areas (compare Figs 1 and 4).
Figure 5 Comparison of average net gain in suitable habitats across species (proportion of area occupied under future change scenario - proportion of area occupied under current conditions) within the study area (SA) and protected areas (PA, including ecological [Ec] network). CL = climate change scenario, CL&HED = scenario of climate change plus high economic development.
DISCUSSION
The habitat suitability modelling suggests a potentially great expansion of the 11 alien invasive plant species in the Ukrainian Carpathians, reaching the subalpine line and occurring along virtually all major linear and/or frequently disturbed habitats in the mountains over the next 100 years. The consistently high permutation importance given to predictors depicting disturbances and the spatial patterns of potential migration into the interior of the UA Carpathians along rivers and roads emphasize the importance of these linear habitats for the establishment of propagating populations of invasives (Prots & Drescher Reference Prots and Drescher2010; Tickner et al. Reference Tickner, Angold, Gurnell and Mountford2011). Although fewer sites are predicted/projected to be at risk of invasion by all nine study species/genera in protected areas (PAs), all PAs are at an increasing risk of being invaded by at least one species. In particular, habitat suitability is projected to increase dramatically at low and medium altitudes by 2050 and 2100, while the ecological network connecting PAs is proposed to cover many of these high-suitability areas.
A crucial step in designing the ecological network was to divide the mountain range into units (approximately 815 × 815 m) of low (0) to high (50) potential to meet conservation targets (Zingstra et al. Reference Zingstra, Seffer, Lasak, Guttova, Baltzer, Bouwma, Walters, Smith, Kitnaes, Predoiu, Prots and Sekulic2009). Potential invasibility is generally concentrated in units of low conservation value (low potential to meet targets), which correspond to areas outside of PAs and the proposed ecological network; these were assigned conservation values of > 10 (Zingstra et al. Reference Zingstra, Seffer, Lasak, Guttova, Baltzer, Bouwma, Walters, Smith, Kitnaes, Predoiu, Prots and Sekulic2009). At the same time, units of highest priority for conservation (value 50) are also relatively more exposed to invasives (Fig. 6).
Figure 6 Invasibility of regions with different conservation values by (a) at least one invasive plant species and (b) all nine species/genera modelled with Maxent under current environmental conditions and future projections. Invasibility is shown as proportion of total units of a specific conservation value (0–50) within the Ukrainian Carpathians in which at least one or all nine invasives can potentially occur (CL&HED = climate change plus high economic development).
Under current conditions, almost 40 % of these units are predicted to contain suitable habitat for at least one study species and over 15 % for all nine species/genera (Fig. 6). These ‘high-conservation-high-potential-invasibility’ units occur particularly along freshwater habitats at low altitudes in the south, south-west and east (Fig. 4; Zingstra et al. Reference Zingstra, Seffer, Lasak, Guttova, Baltzer, Bouwma, Walters, Smith, Kitnaes, Predoiu, Prots and Sekulic2009). This suggests that invasive plants are highly likely to become a permanent feature of ecosystems within the currently proposed network if the observed increasing spread continues (Prots & Drescher Reference Prots and Drescher2010). Because the invasive species might alter the structure and functioning of natural and semi-natural ecosystems (Botta-Dukàt & Balogh Reference Botta-Dukàt and Balogh2008), the potential susceptibility of protected areas and the ecological network to invasion is alarming.
However, the Maxent models used in this study to estimate habitat suitability depict potential establishment sites and not likelihood of invasion, as they do not consider community interactions (such as competitive exclusion) once propagules reach the suitable habitats, or dispersal mechanisms that may facilitate or slow the spread into suitable habitats (Austin Reference Austin2002; Evangelista et al. Reference Evangelista, Sunil, Stohlgren, Jarnevich, Crall, Norman and Barnett2008; Dullinger et al. Reference Dullinger, Kleinbauer, Peterseil, Smolik and Essel2009). In addition, projections assumed static niche preferences of the species, although genetic adaptations and thus potential spread into wider habitats can by no means be ruled out (Lavergne & Molofsky Reference Lavergne and Molofsky2007; Dlugosch & Parker Reference Dlugosch and Parker2008; Clements & Ditommaso Reference Clements and Ditommaso2011).
From potential to probable spread
The issue of leaving out covariates related to dispersal, interactions with native species, and adaptation, all of which in part explain the current distribution of the study species, becomes evident by running a Moran's I correlogram (Fortin et al. Reference Fortin, Dale and ver Hoef2006) on model residuals (1 – predictions on presence locations), which, on average, shows significant spatial autocorrelation (p < 0.001) at distance categories of 1000–10 000 m (1000 m intervals; spatial weight based on inversed distance). Spatial autocorrelation in the model residuals may be explained by missing predictors on adaptation, community interactions or dispersal (Legendre Reference Legendre1993; Dormann et al. Reference Dormann, McPherson, Araújo, Bivand, Bolliger, Carl, Davies, Hirzel, Jetz, Kissling, Kühn, Ohlemüller, Peres-Neto, Reineking, Schröder, Schurr and Wilson2007) and must be investigated in future studies. In general, detailed genetic and demographic studies on the invasive species to approximate future migrations (see for example Genton et al. Reference Genton, Shykoff and Giraud2005; Jongejans et al. Reference Jongejans, Skarpaas and Shea2008) and integrative approaches to estimate the likelihood of establishment given habitat suitability and dispersal/community dynamics (see for example Wadsworth et al. Reference Wadsworth, Collingham, Willis, Huntley and Hulme2000; Nenzén et al. Reference Nenzén, Swab, Keith and Araújo2012; Smith et al. Reference Smith, van Klinken, Seabrook and McAlpine2012) are needed.
Meanwhile, despite a lack of quantitative data on dispersal within the UA Carpathians, the invasive species at the focus of this study have been chosen precisely because they are successful colonizers and long-distance dispersers that use waterways and/or humans for dispersal of propagules and vegetative parts (Williamson et al. Reference Williamson, Pyšek, Jarosik and Prach2005; Botta-Dukàt & Balogh Reference Botta-Dukàt and Balogh2008). In addition, most areas within the study region are relatively easily accessible and connected via an extensive hydrological network and road system that facilitates potential plant invasions. Thus, species-specific traits and area-specific spatial composition indicate that the study species are likely to reach most, if not all, suitable habitats for establishment, including PAs and the ecological network where, despite protection of riparian habitats, natural disturbances (such as floods) may lead to establishment of invasives, given the availability of progagules (Pyšek & Prach Reference Pyšek, Prach, de Waal, Child and Brock1994). Model predictions for current conditions do not extend far beyond areas in which one or several species already occur (compare Fig. 1 and Fig. 2).This supports the hypothesis that the species are highly successful in dispersing across the invaded range.
Similarly, although systematic studies on the genetic composition of populations are needed, many of the species have been chosen for this study due to their high potential for rapid evolutionary adaptations (see Bailey et al. Reference Bailey, Bímová and Mandák2007; Chun et al. Reference Chun, Corre and Bretagnolle2011; Erfmeier et al. Reference Erfmeier, Böhnke and Bruelheide2011). Populations of Heracleum sosnowskyi have increasingly been observed spreading into grasslands and secondary forests away from linear habitats (B. Prots, personal observation). Thus, a precautionary management approach would consider potential spread of the study species and hence the impact on ecosystems beyond current suitable habitats. In addition, future spread of the invasives beyond the areas estimated as suitable by Maxent is likely because the data on predictors used for model calibration stem from a restricted part of the species’ ranges. This leads to an underestimation of the climatic niche (Guisan & Zimmermann Reference Guisan and Zimmermann2000; Franklin Reference Franklin2009) and increases the probability that PAs will be affected more seriously under climate change.
Given the high risk of invasion of suitable habitats, analyses of potential introduction of the invasive species into PAs must be incorporated into strategies to protect biodiversity and must be included in the support of the ecological network (Townsend & Levey Reference Townsend and Levey2005; Zingstra et al. Reference Zingstra, Seffer, Lasak, Guttova, Baltzer, Bouwma, Walters, Smith, Kitnaes, Predoiu, Prots and Sekulic2009).
Future change projections
Although potential range expansions of invasive plant species under climate and land-use change have been shown in several studies (see Ficetola et al. Reference Ficetola, Maiorano, Falcucci, Dendoncker, Boitani, Padoa-Schioppa, Miaud and Thuiller2010; Murray et al. Reference Murray, Stokes and van Klinken2012), projections of potential species distributions must be interpreted with care in this study. First, climate projections used here apply to only one scenario and do not evaluate the range of potential trajectories in regional climate change and thus in potential species distributions. In addition, by adding one value across all climatic pixels for projections, regional differences in climate change, for example between the north-east and south-west, are generalized.
Second, this study does not incorporate detailed scenarios of land-use change and human population development. That is, changes in disturbances/propagule pressure are modelled as intensification (by 10% and 30%) of current spatial patterns for these processes (such as proximity to existing human development). In reality however, changes in disturbances influencing the potential distribution of invasives are expected to be spatially highly dynamic (Verburg et al. Reference Verburg, de Koning, Kok, Veldkamp and Bouma1999; Rounsevell et al. Reference Rounsevell, Reginster, Araujo, Carter, Dendoncker, Ewert, Hause, Kankaanpa, Leemans, Metzger, Schmit, Smith and Tuck2006). Within the study area, increasing illegal forestry practices (Sitko and Troll Reference Sitko and Troll2008; Kuemmerle et al. Reference Kuemmerle, Chaskovskyy, Knorn, Radeloff, Kruhlov, Keeton and Hostert2009), urban and infrastructure development at lower altitudes, and farmland abandonment in the interior of the mountain range (Turnock Reference Turnock2002; Baumann et al. Reference Baumann, Kuemmerle, Elbakidze, Ozdogan, Radeloff, Keuler, Prishchepov, Kruhlov and Hostert2011) are also important processes. While illegal forestry and urban and infrastructure development increase the likelihood of the introduction, establishment and dispersal of invasives (Colautti et al. Reference Colautti, Grigorovich and MacIsaac2006), farmland abandonment may permit the regeneration of natural communities that act as buffers to invasion (Shea & Chesson Reference Shea and Chesson2002). Pauchard et al. (Reference Pauchard, Kueffer, Dietz, Daehler, Alexander, Edwards, Arévalo, Cavieres, Guisan, Haider, Jakobs, McDougall, Millar, Naylor, Parks, Rew and Seipel2009) observed that plant invasions in mountain regions become problematic in open sites. Thus, a road leading from a focal area of current plant invasion to a deforested patch may provide a corridor for invasives to the open area. In PAs, projected increasing investments into tourism (Webster et al. Reference Webster, Holt and Avis2001) may create novel suitable habitats for invasives beyond current linear trajectories (see for example Dickens et al. Reference Dickens, Gerhardt and Collinge2005).
CONCLUSION
Regional climatic and land-use projection models at a finer resolution are urgently needed for the UA Carpathians in order to develop more detailed and realistic species distribution projections. Nevertheless, our simple projections clearly illustrate that, given the apparent dispersal success of the 11 study species, invasion of areas of high conservation value is likely if current trends continue, and eradication of invasive species is highly unlikely. Current conservation planning in the UA Carpathians must be amended to include the long-term presence and management of invasive species in PAs and the ecological network.
ACKNOWLEDGEMENTS
We thank Dr Ivan Kruhlov (Lviv), Dr Tobias Kuemmerle (Berlin), Dr Patrick Hostert (Berlin), Dr Yuriy Andreychuk (Lviv), and Dr Taras Yamelynets (Lviv) for allowing us to use their grid and vector maps on the Ukrainian Carpathians and for providing advice during our studies.
