Modelling plant invasion pathways in protected areas under climate change: implication for invasion management

Global climate change may enable invasive plant species (IPS) to invade protected areas (PAs), but plant invasion on a global scale has not yet been explicitly addressed. Here, we mapped the potential invasion pathways for IPS in PAs across the globe and explored potential factors determining the pathways of plant invasion under climate change. We used species distribution modelling to estimate the suitable habitats of 386 IPS and applied a corridor analysis to compute the potential pathways of IPS in PAs under climate change. Subsequently, we analysed the potential factors affecting the pathways in PAs. According to our results, the main potential pathways of IPS in PAs are in Europe, eastern Australia, New Zealand, southern Africa, and eastern regions of South America and are strongly influenced by changes in temperature and precipitation. Protected areas can play an important role in preventing and controlling the spread of IPS under climate change. This is due to the fact that measures are taken to monitor climate change in detail, to provide effective management near or inside PAs, and to control the introduction of IPS with a high capacity for natural dispersal. A review of conservation policies in PAs is urgently needed.


Introduction
Invasive plant species (IPS) have the potential to threaten global biodiversity in the face of global climatic changes (Kalusová et al., 2013).However, the relationship between IPS and climate change is complex (Hellmann et al., 2008).Different patterns of climate change may result in different distributions of species on regional scales (Bradley, 2010), and climatic factors are among the primary factors determining the overall distribution patterns of IPS due to potential synergic effects (Bradley, 2010;Bai et al., 2013).Such effects are related to the influence of climate factors on the global terrestrial net primary production, water availability, plant diseases, and the large-scale sexual reproduction of IPS (Rosenzweig et al., 2001;Hedhly et al., 2009;Kaser et al., 2010;Zhao and Running, 2010).Climate change can promote the expansion of IPS into non-native ranges, decrease ecosystem stability, and threaten native plant diversity (Hellmann et al., 2008;Richardson and Rejmánek, 2011;Bai et al., 2013).The increasing expansion of IPS can facilitate their establishment and naturalization (Thuiller et al., 2005;Wilson et al., 2009;Kalusová et al., 2013;Donaldson et al., 2014).Human activities, such as agriculture, forestry, and horticulture, also increasingly interact with climate change and affect the invasion pathways of IPS (Reichard and White, 2001;Donaldson et al., 2014).Intentionally introducing IPS into new areas may play an increasingly important role in the invasion of introduced plant species, especially due to the potential synergic effects of climatic change and human activity.
Niche conservatism in space and time is one of the key assumptions for analysing climate change impacts on IPS movements (Petitpierre et al., 2012).Petitpierre et al. (2012) showed that IPS can grow and survive in regions with environmental conditions similar to those in their native ranges.However, some studies indicated that the climatic niches of IPS may shift between native and invasive ranges and that IPS have the ability to occupy new environmental niches that differ from their native ranges (Early and Sax, 2014;Dellinger et al., 2016;Wan et al., 2016a, b;Wang et al., 2017).This implies that IPS expand and establish under the same environmental conditions in native and invasive ranges in the above-mentioned situations (Petitpierre et al., 2012;Early and Sax, 2014;Wan et al., 2016a, b;Wang et al., 2017).Hence, IPS can invade ecologically fragile protected areas (PAs) with environmental conditions similar to those of the native and invaded environments of these species (Foxcroft et al., 2013;Vicente et al., 2013;Hulme et al., 2014).Protected areas often contain diverse fauna and flora and unusual geological features; they therefore play a crucial role in biodiversity conservation (Foxcroft et al., 2013).Previous studies have shown that IPS have the potential to invade regional PAs under the scenario of a changing climate (Foxcroft et al., 2013), consequently altering species composition, changing the ecosystem functions, and even causing PAs to lose their function in the conservation of endangered species (Pauchard and Alaback, 2004;Hulme et al., 2014).Foxcroft et al. (2011) have suggested that PA boundaries could limit the spread of IPS into the interior of the PAs.However, IPS can pass through the boundaries of PAs and invade PAs via human activities (e.g.transportation, trade and forestry) and natural dispersion (McConnachie et al., 2012;Foxcroft et al., 2013;Donaldson et al., 2014).The connectivity of PAs, e.g. through corridors, can even offer additional routes for the spread of IPS (Vicente et al., 2013;Donaldson et al., 2014).Therefore, it is essential to evaluate whether climate change affects the movement of IPS through PAs and whether PAs can act as barriers to the invasion of non-native species.
Species distribution modelling (SDM) is widely used to predict suitable habitats of plant species based on niche conservatism on the global scale; it can also be used to analyse the complex relationships between IPS and PAs (Václavík and Meentemeyer, 2009;Petitpierre et al., 2012;Vicente et al., 2013;Kolanowska, 2014;McConnachie et al., 2015).Species distribution modelling coupled with geographic information systems (GIS) can be applied to generate spatially explicit models of the potential of IPS to invade PAs and of the invasion pathways under climate change (Thuiller et al., 2005;Vicente et al., 2013;Donaldson et al., 2014).Recent studies have shown that IPS can invade PAs under climate change at the regional scale (Vicente et al., 2013;Thalmann et al., 2015).However, our current understanding of the relationship between climate change, PAs, and invasion pathways of IPS on a global scale is extremely limited.We hypothesized that PAs play an important role in the prevention and control of IPS under climate change.To test this hypothesis, we used two indices: (1) the potential factors of IPS that invade PAs under climate change and (2) the overlap of potential pathways of IPS and PAs, including pathways due to the intentional introduction of IPS.
We selected 386 IPS from the Invasive Species Specialist Group (ISSG) list and used Maxent to model the potential distribution of these 386 IPS on a global scale (http: //www.issg.org/database/species/List.asp).We then used corridor analysis in GIS to model the invasion pathways of IPS affected by climate change and compared the overlap of potential pathways of IPS and global PAs (Thuiller et al., 2005).Finally, we analysed the potential factors affecting IPS pathways in PAs.

Protected area data
The International Union for Conservation of Nature (IUCN) has developed management categories based on the management objectives of PAs (Dudley, 2008).These categories provide international criteria for defining PAs and encourage conservation planning according to the management requirements of PAs on a global scale.We downloaded a global map of the IUCN I-VI PAs from the World Database on Protected Areas (WDPA; https://www.protectedplanet.net/).For our analysis, we selected over 23 852 PAs with areas over 10.0 × 10.0 km 2 , as they have a large enough ecological landscape to meet the requirements for studying the potential extent of pathways of IPS and to match the relatively coarse resolution of bioclimatic data.

Bioclimatic data
Current and future data used for modelling were 5.0 arcmin for the environmental layer input of the SDM.Nineteen common bioclimatic variables were downloaded from the WorldClim database (averages from 1950-2000 were used as current bioclimatic variables; detailed information in www.worldclim.org;Bellard et al., 2014).Bioclimatic variables with Pearson's correlation coefficients between 0.7 and −0.7 (using Spearman non-normal statistic) were removed to eliminate the negative effect of multicollinearity on SDM adjustment.The remaining four bioclimatic variables (Table S1 in the Supplement) influence the distribution and physiological performance of the IPS (Gallagher et al., 2013).These four bioclimatic variables imply the mean and standard deviation of current and future climates (Hijmans and Graham, 2006).We relied on future climate data from the Intergovernmental Panel on Climate Change (IPCC), Fifth Assessment Report (AR5) as a reference for modelling the changing trends of IPS invasion (http://www.ipcc.ch/).To model the potential future distribution of IPS in the 2080s (2071-2099), we used an average map of four global climate models (GCMs; i.e. bcc_csm1_1, csiro_mk3_6_0, gfdl.cm3, and mohc_hadgem2_es) and two greenhouse gas emission scenarios, the Representative Concentration Pathways (RCPs) 4.5 (mean: 780 ppm; range: 595 to 1005 by 2100) and 8.5 (mean: 1685 ppm; range: 1415 to 1910 by 2100), representing the low-and high-gas-concentration scenarios, respectively (http://www.ccafs-climate.org/).

Species data
We used 386 IPS with more than 50 records from the Invasive Species Specialist Group (ISSG) list of IUCN as a representative set of the IPS around the world (Wisz et al., 2008;Merow et al., 2013;García-Roselló et al., 2014;Duputié et al., 2014;Table S2).Occurrence data, in particular geographic coordinates, for each IPS were obtained from the Global Biodiversity Information Facility (GBIF; www.gbif.org),which is one of the largest online providers of distribution records.We removed duplicate occurrences of recorded data for species in 5.0 arcmin grid cells (10 km at the equator) to avoid any georeferencing errors (Merow et al., 2013).All occurrence localities were checked using Google Earth and ArcGIS 10.2 (Esri; Redlands, CA, USA) to determine whether IPS were distributed in reasonable ranges based on the ISSG; obvious errors were removed (Bellard et al., 2014).

Species distribution modelling
We used a machine-learning technique called "maximum entropy modelling" (hereafter referred to as "Maxent") to model suitable habitats for each IPS in current lowand high-gas-concentration scenarios, based on current presence-only species records (Phillips et al., 2006;Smith et al., 2012;Merow et al., 2013;Kolanowska, 2014).We ran Maxent using Maxent software version 3.4.0downloaded from http://biodiversityinformatics.amnh.org/open_source/maxent/.Maxent can express a probability distribution where each grid cell has a predicted suitability of conditions for the species (http://biodiversityinformatics.amnh.org/open_source/maxent/).The set of all grid cells was regarded as the possible distribution space of maximum entropy (namely, that is most spread out, or closest to uniform; Phillips et al., 2006).For the map cells predicted using Maxent, cells with values of 1 had the highest degree of habitat suitability, while cells with values of 0 had the lowest (Elith et al., 2011).Habitat suitability was determined based on the climatic similarity to sites where the species already occur (Papeş and Gaubert, 2007;Holcombe et al., 2007).Our modelling sets were the same as described in Wan et al. (2016a).
The predictive precision of Maxent was based on the area under the curve (AUC) of the receiver operation characteristic (ROC), which regards each value of the prediction result as a possible threshold; the corresponding sensitivity and specificity were obtained through calculations.The AUC ranges from 0.5 (lowest predictive ability or not different from a randomly selected predictive distribution) to 1 (highest predictive ability).Models of each species with values above 0.7 were considered useful in our study.Only one species, Poa pratensis, had an AUC value below 0.7 and was hence not considered in downstream analyses (Hijmans, 2012; Table S2).The AUC values of the other 385 species were above 0.7, indicating adequate model performance.Fi-nally, we combined the Maxent results of the 385 IPS with the current and future climatic data to develop maps of suitable current and future IPS habitats.

Modelling potential invasion pathways
First, IPS distribution pathways were estimated using a corridor analysis with the input variables set to two resistance layers (ESRI, 2014) to identify the resistance of passing grids.Resistance layers indicated the movement cost of IPS across spatial and temporal scales.We considered the grids with low values of suitable habitats as resistance grids and the grids with high values of suitable IPS habitats as non-resistance grids.The purpose of a corridor analysis (the corridor function in ArcGIS 10.2; Esri; Redlands, CA, USA) is to build a new pathway of minimum resistance of passing grids; this method was implemented between two resistance layers, using the corridor function in ArcGIS 10.2 (Esri; Redlands, CA, USA; Rabinowitz and Zeller, 2010).Suitable habitats of low and high concentration scenarios were considered as potential pathways for IPS distribution under current and future conditions (Morato et al., 2014).This dispersal scenario simulation of the corridor analysis served to identify connected or unconnected pairs of patches, namely invasion locations from the current to the future (including low-and high-gasconcentration scenarios; Morato et al., 2014).
Subsequently, we identified PAs which covered large areas of potential movement pathways of IPS.We used the following equation to compute the potential pathways for IPS in PAs (Alagador et al., 2011;Bellard et al., 2014): where C j represents the ability of PAs to support the potential pathway of IPS, n is the number of species in PA j , and P j,n is the average habitat suitability of cells for species n in PA j .The larger values represented the stronger ability of PAs to support the potential pathways of the IPS (i.e.index of potential pathways in PAs).A linear-regression analysis was used to compute the ability of the PAs to support the potential pathways of the IPS under low-and high-gas-concentration scenarios (Peterson et al., 2008).

Potential factors determining invasion pathways in protected areas
We compared the ranges of potential pathways of IPS in PAs based on the Jenks optimization method, a clustering method designed to determine the best arrangement of values into different classes, based on the reduction and maximization of variance within and among classes.The purpose was to determine the impact of climatic change on potential IPS pathways.The most important climatic variables (contribution percentage to Maxent modelling larger than 15 %) were www.web-ecol.net/17/69/2017/Web Ecol., 17, 69-77, 2017 determined using a Jackknife test in the programme Maxent (Merow et al., 2013).We computed the average variances of the most important climatic variables that are most closely related to habitat suitability for IPS in each PA and classified the variance into three groups (low, medium, and high degrees of variance of climate change).Subsequently, we applied a two-way ANOVA to compare the differences in the potential IPS pathways between the three groups, based on the most important temperature and precipitation variables (Baz et al., 2009).All analyses were conducted in JMP 11.0 (SAS Institute).

Potential invasion pathways in protected areas
Our results show that IPS are most likely to spread across Europe, eastern Australia, New Zealand, southern Africa, and the eastern regions of South America in both low-and high-gas-concentration scenarios (Fig. 1).In the PAs, there was a significant relationship of the potential IPS movement pathways between low-and high-gas-concentration scenarios (P < 0.001; Fig. 2).Hence, we used the PA map of potential IPS pathways in the low-gas-concentration scenario as the representative map for the potential IPS pathways in PAs.
The degree of the potential pathways of IPS, including maximum, minimum, and square deviation values in PAs, was similar to that of the global potential pathways (Table 1).The maximum degree of potential IPS pathways was slightly lower in PAs than across the globe (Table 1).Minimum and average values of potential pathways of IPS were higher in PAs than in a global context in both low-and highconcentration scenarios (Table 1).

Potential factors determining invasion pathways in protected areas
The most important climatic variables that might facilitate the spread of IPS were temperature changes, including annual mean temperature (Bio1), temperature seasonality (Bio4), and annual precipitation (Bio 12; Table S2).Annual mean temperature coupled with annual precipitation had synergic effects on potential invasion pathways in PAs in a lowgas-concentration scenario (two-way ANOVA test between the three groups: P < 0.05; Table 2).With increasing variance in annual mean temperature, the potential movement of IPS decreased in both low and high-concentration scenarios (ANOVA test between the three groups: P < 0.05; Table 2; Fig. 3a).Also, the variation of annual precipitation may have significant negative effects on the potential movement of IPS in both low-and high-gas-concentration scenarios, except for the medium vs. high variance of annual precipitation in the high-gas-concentration scenarios (ANOVA test between the three groups: P < 0.05; Table 2; Fig. 3b).

Potential factors determining invasion pathways in protected areas
Our results suggest that the changes in annual temperature and precipitation determine the potential pathways for invasive plant species (IPS) inside and outside of protected areas (PAs).For PAs with numerous potential pathways for IPS, suitable habitats seem to be crucial for colonization.Our results suggest that the extent of potential pathways for IPS in PAs is likely to increase with temperature.The large variability in precipitation in the PAs may lead to numerous potential pathways.The results of our study also indicate that, in order to be able to spread, IPS require appropriate temperatures with pronounced precipitation.Some other studies have also shown that climate change can lead to the creation of suitable habitats for the invasion of IPS in PAs (Pickering et al., 2011;Vicente et al., 2013;Thalmann et al., 2015).The connectivity of these suitable habitats can then provide important corridors for IPS movements (Vicente et al., 2013).Moreover, IPS often have the ability to spread rapidly in PAs (Foxcroft et al., 2013), and extreme weather events, such as large seasonal differences in temperature and precipitation within a year, can facilitate the formation of pathways for IPS in PAs (Bradley et al., 2010;Diez et al., 2012).Furthermore, the proceeding globalization facilitates the spread of IPS, as international commerce develops new trade routes, markets, and products (Perrings et al., 2005;Bradley et al., 2010;Seebens et al., 2015).For example, emerging economies in mega-diverse regions expect a strong increase in the number of naturalized plants by 2030 (Seebens et al., 2015).Anthropogenic activities, such as agriculture and forestry in regions surrounding PAs, further promote the introduction of IPS (Donaldson et al., 2014;Melin et al., 2014;McConnachie et al., 2015).Intensive anthropogenic activities inside PAs and in adjacent areas greatly increase the number of pathways for IPS (Spear et al., 2013;Mc-Connachie et al., 2015).As habitat suitability for IPS shifts with climate change, the natural dispersal of IPS, as opposed to the anthropogenic dispersal, could also promote invasion into PAs (www.issg.org).According to Colautti and Bar-  SDM has some uncertainties in the distribution predictions.We need to take potential niche shifts of IPS between native and invasive ranges into consideration and use more occurrence records covering the niches of native and invasive ranges to model the distributions of IPS under climate change (Early and Sax, 2014;Dellinger et al., 2016;Wan et al., 2016a, b;Wang et al., 2017).With accelerating economic globalization and rapid climate change, a risk analysis of IPS on regional and global scales and intensive studies of the dispersal methods of IPS are crucial.

Preventing and controlling IPS in protected areas
There was a high degree of overlap between IPS pathways for PAs and in a global context, particularly in Europe, eastern Australia, New Zealand, southern Africa, and eastern regions of South America, and the degrees of pathways in PAs were close to those of the global potential pathways of IPS.Previous studies have shown that potential pathways of IPS may only become apparent in the late stages of invasion (Meier et al., 2014;Donaldson et al., 2014), and our findings indicate that the potential pathways of IPS in PAs play an important role in the intentional introduction of IPS worldwide.Thus, IPS may not only impact biodiversity and habitat quality within PAs (Foxcroft et al., 2013;Hulme et al., 2014), but also have the potential to easily spread to adjacent PAs (Thuiller et al., 2005;Kalusová et al., 2013).Hence, PAs may promote the development of potential pathways for IPS and then increase the risk of invasion on a global basis, given that PA management is not being enhanced.Preventing and con-   trolling pathways for IPS and thus reducing their substantial impacts are therefore a crucial component of conservation management in PAs (Van Wilgen et al., 2012).
Measures should be taken to prevent the invasion of PAs in order to maximize their capacity to withstand colonization by IPS (Foxcroft et al., 2011).Preventing IPS from colonizing PAs can decrease the population sizes of IPS (Foxcroft et al., 2013).The challenge for environmental managers will be to minimize the opportunities for IPS to be introduced into new areas, especially in the face of climate change.Adequate management has the potential to control IPS movements and therefore to prevent invasion into new areas (Hellmann et al., 2008;Lockwood et al., 2012).Such prevention depends on early detection and mechanisms to avoid introductions of IPS (Meier et al., 2014;Donaldson et al., 2014), whereas IPS control is related to measures to remove or eradicate invasive species in a particular area (Foxcroft et al., 2011(Foxcroft et al., , 2013;;Meier et al., 2014).Hence, a review of the conservation policies of PAs is of utmost importance (Van Wilgen et al., 2012;Cuddington et al., 2013).For example, European legislation plays a crucial role in developing policies to protect the environment and meeting its objectives for sustainable development (http://jncc.defra.gov.uk/page-1372).Europe is the key region of prevention and control of plant invasions under climate change.We need to integrate climate change and the pathways of plant invasion into the conservation policies for the management of PAs.It is therefore necessary to study local dispersal methods of IPS, to control the natural dispersal of IPS near or inside PAs, and to prevent intentional introduction of IPS (Van Wilgen et al., 2012).

Conclusions
We used a novel approach to model potential invasion pathways of IPS under climate change and explored the spatial overlap between plant invasion pathways and PAs.Global changes in temperature are likely to control potential pathways of IPS in the future, and the high degree of overlap between PAs and these potential global invasion pathways might pose a major problem for biological conservation.Europe, eastern Australia, New Zealand, southern Africa, and eastern regions of South America are the key regions of IPS prevention and control.The changes in annual temperature and precipitation are likely to promote the movement of IPS in the future.Our results thus suggest that measures blocking pathways for IPS in PAs could effectively control the spread of IPS worldwide.Further research is needed to complement our study, such as modelling plant invasions with greater spatial resolution (Cuddington et al., 2013;Hulme et al., 2014).
Data availability.Data are available on request.

Figure 1 .
Figure 1.Maps of the global potential pathways of invasive plant species (IPS) in protected areas (PAs) under climate change.(a) Global potential pathways of IPS in a low-gas-concentration scenario; (b) Global potential pathways of IPS in a high-gas-concentration scenario.Black arrows represent the highlighted potential pathways of IPS in PAs under climate change.

Figure 2 .
Figure 2. Relationship of the pathway degrees for invasive plant species (IPS) in protected areas (PAs) between low-and high-gasconcentration scenarios.

Figure 3 .
Figure 3. Ranges of potential pathways for invasive plant species (IPS) in protected areas under climate change.CC is the degree of variation of climate change.CS is the gas concentration scenario.Bars represent standard error; vertical axes represent the extent of potential pathway degrees for IPS in protected areas.

Table 1 .
Comparison of pathways for invasive plant species (IPS) between protected areas (PAs) and global regions based on the map cells.Min is the minimum pathway values for IPS between PAs and global regions, max is the maximum, mean is the mean average, and SD is standard deviation.

Table 2 .
Effects of annual mean temperature and annual precipitation on potential invasion pathways in protected areas (PAs) in both lowand high-gas-concentration scenarios based on two-way ANOVA.
Bio1 is annual mean temperature, Bio12 is annual precipitation, CS is the gas concentration scenario.F , P and degree of freedom are based on a two-way ANOVA.