Climate change amplifies plant invasion hotspots in Nepal

Climate change has increased the risk of biological invasions, particularly by increasing the climatically suitable regions for invasive alien species. The distribution of many native and invasive species has been predicted to change under future climate. We performed species distribution modelling of invasive alien plants (IAPs) to identify hotspots under current and future climate scenarios in Nepal, a country ranked among the most vulnerable countries to biological invasions and climate change in the world.

Invasive alien plants (IAPs; sensu Pyšek et al., 2014) pose the greatest threats to natural ecosystems, human health, economy, agriculture and fisheries (Pimentel, Zuniga, & Morrison, 2009;Vilà et al., 2010;Vilà & Hulme, 2017). The threat and loss constituted by invasive species are exacerbated by climate change through multiple mechanisms including the removal of climate barriers for establishment and the spread of many invasive species (Bradley, Blumenthal, Wilcove, & Ziska, 2010;Hellmann, Byers, Bierwagen, & Dukes, 2008). For example, Petitpierre et al. (2017) have shown that the upslope spread of the lowland IAPs to mountains is limited by low temperatures in the mountains. Yet, climate change will shift regions of optimal suitability for these IAPs from lowlands to highlands.
Invasive species have a greater capacity to shift their niches more rapidly than native species, and they are more likely to adapt to new climatic conditions faster (Dukes & Mooney, 2004;Hellmann et al., 2008). The climatic niche shift has been demonstrated as one of the mechanisms used by IAPs to spread rapidly into introduced ranges (e.g., Treier et al., 2009;Gallagher, Beaumont, Hughes, & Leishman, 2017). Furthermore, the IAPs benefit from global warming and atmospheric CO 2 enrichment more than native plants (Liu et al., 2013;Verlinden & Nijs, 2010). Therefore, an integrated understanding of biological invasions and climate change is necessary for the management of IAPs.
Projected changes in climate in the future may influence the distribution of many native and invasive species (Bellard, Bertelsmeier, Leadley, Thuiller, & Courchamp, 2012;Walther et al., 2009). At the species level, climate change causes range expansions of many invasive species and contraction to few (Bellard et al., 2013;Bradley et al., 2010). Therefore, these potential changes in distribution need to be incorporated into the management and conservation of ecosystems and biodiversity in the face of biological invasions and climate change (O'Donnell et al., 2012). A better understanding of threats and an ability to accurately predict the impacts of climate change on species distribution are necessary to make an informed decision for biodiversity conservation (Pimm et al., 2005). This will help to minimize the threat of invasive species into the future and support effective conservation efforts. Although the impact of climate change on the distribution of multiple invasive species is known in developed regions such as Australia, North America and Europe (Allen & Bradley, 2016;O'Donnell et al., 2012), little is known about how the distribution of invasive species will change with future climatic changes in developing countries such as Nepal.
With the estimated annual cost of US$ 1.4 billion, due to biological invasions to Nepal's agriculture sector, Nepal is ranked among the topmost countries (ranked 3 rd out of 124 countries) in terms of invasion threats to agriculture sectors (Paini et al., 2016).
The concentration of vertebrate species threatened by biological invasions is also high in the Indian subcontinent including Nepal (Bellard, Genovesi, & Jeschke, 2016). Currently, there are 241 alien plants and animals in Nepal and 45 of them are considered invasive (Shrestha, Budha, Wong, & Pagad, 2018). These invasive species can be found from lowland plains in the south to hills and mountains in the north. Globally, mountain ecosystems are generally less invaded compared to the surrounding lowlands (McDougall et al., 2011). However, the intensity of biological invasions is likely to increase in future with changing climate and increasing anthropogenic disturbances in the mountains (Pauchard et al., 2016;Petitpierre et al., 2017). In Nepal, most of the IAPs are found below 2,000 m in elevation (Shrestha, 2016) but recent studies based on field observations and models suggest that some of these IAPs are already expanding their ranges into new geographic locations at a higher elevation (Lamsal, Kumar, Aryal, & Atreya, 2012;Thapa, Chitale, Rijal, Bisht, & Shrestha, 2018).
The severity of threats to Nepal's economy and ecosystems from biological invasions is considered in national conservation policies and sectoral conservation strategies such as the Plant Protection Act (2007), National Biodiversity Strategy and Action Plan (2014), Forestry Sector Strategy (2016-2015 and National Ramsar Strategy and Action Plan (2018)(2019)(2020)(2021)(2022)(2023)(2024). However, the implementation of these policies and strategies is very poor, partially because of the lack of scientific knowledge required to control invasive alien species (MFSC, 2014;Shrestha et al., 2015). National Biodiversity Strategy and Action Plan (2014-2020) has identified priorities of actions for the management of invasive alien species that includes, among others, research and prioritization of problematic IAPs (MFSC, 2014).
To this end, this study is an important contribution to enhance the knowledge and understanding of invasive species by identifying potentially suitable niches for 24 IAPs (out of reported 26 species) in Nepal under current and future climate using an ensemble of species distribution models.
We also performed a hotspot analysis to identify areas suitable for a maximum number of IAPs. In biodiversity conservation, the concept of a biodiversity hotspot-an area with high species richness, endemism and threatened taxa-is well established (e.g., Myers, Mittermeier, Mittermeier, Fonseca, & Kent, 2003). Hotspot analysis provides a framework for cost-effective conservation programmes, thus helping to prioritize conservation efforts (Myers, 2004). Similarly, hotspot analysis can help streamline management efforts in a way to prevent, eradicate and control the maximum number of invasive species at the lowest cost possible (Adhikari, Tiwary, & Barik, 2015;O'Donnell et al., 2012). Attempts at applying the hotspot concept to biological invasions were made at a global scale (Drake & Lodge, 2013), according to country (Australia by O'Donnell et al., 2012, United States by Allen & Bradley, 2016, India by Adhikari et al., 2015 and according to the local scale (Corangamite Catchment in Australia by Catford, Vesk, White, & Wintle, 2011), using different tools and techniques. In this study, we used the invasion hotspot approach to identify regions with high concentrations of potentially suitable niches for multiple IAPs under current and future climate in Nepal. We also examined changes in invasion hotspots according to different land covers, climatic zones, physiographic regions, ecoregions and federal states. This study provides the first-ever comprehensive national-level assessment of biological invasions using occurrence data of most of the IAPs found in Nepal. By highlighting the geographic hotspots for IAPs, our results provide wide-ranging evidence of current and future risks constituted by IAPs at the national scale. The results of this study will be important when considering cost-effective strategies for managing IAPs and will support long-term biodiversity conservation and sustainable development goals in Nepal.

| Invasive alien plants' occurrence data
We selected 24 out of the total 26 species of IAPs reported in Nepal (Shrestha, 2016;Tiwari, Adhikari, Siwakoti, & Subedi, 2005) for the modelling exercise ( Figure 1). Location information sufficient for the distribution modelling of the remaining two species (Myriophyllum aquaticum and Spergula arvensis) is not yet available.
The description (name, family, native origin, functional group, distribution range and mode of dispersal) of the selected species is given in Table 1. Four of the modelled species (Chromolaena odorata, Eichhornia crassipes, Lantana camara and Mikania micrantha) are among the 100 of the world's worst invasive species (Lowe, Browne, Boudjelas, & DePoorter, 2017). Sixteen of the modelled species are considered highly problematic species by the local people of Nepal due to their negative impacts on agriculture, local livelihood and natural ecosystems . Modelled species were introduced either deliberately for ornamental purpose (e.g., L. camara, E. crassipes) or accidentally (e.g., Ageratina adenophora, Parthenium hysterophorus) to Nepal at various times.
The occurrence data were collected through field surveys by experts in various localities of Nepal at different times from 2013 to 2018 (Shrestha, 2014;Shrestha, Joshi et al., 2018;Shrestha, Kokh, & Karki, 2016;Siwakoti et al., 2016). The number of occurrence locations for each species ranged from 25 (Leersia hexandra) to 1,910 (Bidens pilosa). Species distributional data often display spatial autocorrelation which has implications for predicting species occurrences under changing environmental conditions (Dormann, Grime, & Thompson, 2000;Dormann, 2007). We removed the multiple presence locations in the same grid of ~1 km 2 spatial resolution and retained only one unique record per grid by applying spatial filtering F I G U R E 1 Occurrence of 24 invasive alien plants in Nepal. Each dot represents geographic coordinates of the species using sdmtoolbox 2.3 (Brown, 2014). We also checked differences in Moran's I index values after removing the multiple records using the ape package in r. This approach reduces spatial autocorrelation, which could lead to overfitting the models and therefore reduce model performance (Boria, Olson, Goodman, & Anderson, 2014).
After spatial filtering, the total number of occurrence locations comprising 24 species was reduced from 17,682 records to 10,951, and these were then used for modelling.
Bioclimatic variables with high correlation (Pearson's correlation coefficients r > 0.70) were removed to reduce multicollinearity (Rogerson, 2001). The VIF values of the resulting predictor variables were less than 5. A VIF greater than 10 signals a collinearity problem (Chatterjee & Hadi, 2015). The remaining seven bioclimatic variables, namely annual mean temperature, mean diurnal range [mean of monthly (max tempmin temp)], isothermality, temperature annual range (max temperature of warmest month-min temperature of coldest month), precipitation in the driest month, precipitation in the warmest quarter and precipitation in the coldest quarter, were used as predictors to model the current distribution of the selected IAPs.
We also modelled distributions of the 24 IAPs for future climatic conditions. We used projected bioclimatic variables for the period  Table S1.

| Species distribution modelling
Species distribution modelling is an approach that predicts the distribution of a species across geographic space and time using the correlation between the geographic occurrence or abundance of a species and corresponding environmental conditions (Elith & Leathwick, 2010). This approach has been used in studies of biogeography, conservation biology, ecology, palaeoecology and wildlife management for more than a decade (Araújo & Guisan, 2006) and forecasts the range shifts of species under future climate change scenarios (Beaumont, Pitman, Poulsen, & Hughes, 2007;Wiens, Stralberg, Jongsomjit, Howell, & Snyder, 2009) including invasive species (Bellard et al., 2013). Various species distribution modelling tools such as statistical regression, machine learning and geographic extrapolation are in current use to model species distribution (Elith et al., 1999). The performance of various algorithms available for species distribution modelling varies significantly (Elith, Kearney, & Phillips, 2006). An ensemble modelling of species distributions involves simulations across more than one set of initial conditions, model classes, model parameters and boundary conditions (Araújo & New, 2007). BIOMOD (biomod2 package in r) is a platform for ensem- As these models required background data (e.g., pseudo-absence) and the actual absence data were unavailable, we used 10,000 pseudo-absences selected randomly outside a buffer of 10 km from the presence points by following Barbet-Massin, Jiguet, Albert, and Thuiller (2012). The models were calibrated by using 70% of the occurrence points (presence and pseudo-absence) as training data and evaluated by using the remaining 30% as testing data (Araújo, Pearson, Thuiller, & Erhard, 2005). We repeated the process of pseudo-absence generation three times and three evaluation runs per species, resulting in a total of 72 models per species (eight models, three evaluation runs and three pseudo-absence selection procedures) under each climate scenario.
From the 72 models per species, we built ensemble models using a weighted-mean approach in which weights are awarded for each model proportionally to their evaluation metrics scores; hence, the discrimination is fair in this approach (Marmion, Parviainen, Luoto, Heikkinen, & Thuiller, 2000). Only the models with good predictive accuracy (TSS > 0.6 and AUC > 0.8) were used to build an ensemble from the projection outputs (Bellard et al., 2013;Gallien, Douzet, Pratte, Zimmermann, & Thuiller, 2010;Thuiller et al., 2009). Binary maps (suitable and unsuitable) were produced using the optimal threshold that maximizes the TSS score as a cut-off value, which then converted the projected occurrence probabilities during the cross-validation procedure (Allouche et al., 2006;Liu, White, & Newell, 2016;Marmion et al., 2000). This threshold is unaffected by the prevalence of species occurrence and favours sensitivity

| Invasion hotspot map
We conducted a hotspot analysis ( We calculated changes in the areas of so-called invasion hotspots under current and future climate with respect to climatic zones, land cover, ecoregions, physiographic regions and federal states. We used publicly available maps of ecoregions (Olson et al., 2012), land cover (Uddin et al., 2015), physiography and administrative units. A layer of climatic zones was created by using digital elevation model as tropical (<1,000 m asl), subtropical (1,000-2,000 m asl), temperate (2,000-3,000 m asl), subalpine (3,000-4,000 m asl) and alpine (>4,000 m asl) following Shrestha (2008).

| RE SULTS
The model performance was evaluated by the scores of two (AUC and TSS) performance matrices ( Figure S2). Climate change will increase both the extent and the intensity (invasion hotspots) of the climatically suitable regions for IAPs in Nepal ( Figure 2). Under the current climate, around 59,700 km 2 (40%) of the country were predicted as suitable for IAPs, while 33,600 km 2 (23%) were identified as invasion hotspots in Nepal. The niche extent and invasion hotspots will expand by 2% and 5%, respectively, suggesting an increase in the potential niche of IAPs from current to future climate.
Changes in the invasion hotspots were evident in different climatic zones, ecoregions, land covers, physiographic regions and federal states in Nepal (Figure 3). The maximum increase in the area of invasion hotspots was observed in the tropical zone, which is situated below 1,000 m (+2,747 pixels, 8%). Under the current climate, there were no invasion hotspots in the subalpine region (3,000-4,000 m asl); however, the hotspots will expand towards temperate and subalpine regions in the future, indicating an expansion towards higher elevation regions under future climate. Out of the ten ecoregions, eight coincided with the invasion hotspots.

| D ISCUSS I ON
To our knowledge, this study of modelling the distribution of 24 IAPs presents the most comprehensive analysis of biological invasions in Nepal. To date, most research either has focused on a handful of species  or has been limited to smaller geographic area of Nepal . We identified the geographic areas with different land cover, ecoregions, physiography and climatic zones that are climatically suitable for invasions using a novel approach of invasion hotspots. As consistent with the previous studies Thapa et al., 2018), our results show that changing climate will create additional climatically The climatic condition of central Nepal, characterized by subtropical climate that is under cropping and forests, has the maximum suitable areas for a majority of IAPs under current climate. As the climate changes, a new suit of habitats will emerge that may be suitable for alien species (Hellmann et al., 2008). Climate change facilitates dispersal, introduction and naturalization of alien species as well as reduces the resilience of local ecosystems to alien species (Walther et al., 2009). Our models also predict that additional suitable areas for IAPs are expected to emerge in the higher elevation TA B L E 2 Change in the climatically suitable niches (km 2 ) of the invasive alien plant species  showed that IAPs are currently invading higher elevation regions as compared to the past (Shrestha, Shabbir, & Adkins, 2015;Tiwari et al., 2005). At high-elevation regions, the impact of climate change is likely to be more severe than in low elevation, as the magnitude of temperature change is greater in those areas (Shrestha, Gautam, & Bawa, 2012). Our analysis showed that the maximum increase in the area of climatically suitable niches occurred at lower elevations F I G U R E 3 Change in the invasion hotspots for 24 invasive alien plants of Nepal between current and future climate (for 2050 under RCP 6.0). (a) Elevation bands, (b) land use types, (c) ecoregions, (d) physiographic regions and (e) federal states (Nepal has been recently divided into seven federal states and numbered from 1 (east) to 7 (west). Official names of the states have yet to be declared) (below 2,000 m) while the proportion of change in suitable niches is greater at higher elevations (above 2,000 m). The creation of climatically suitable regions for IAPs in the high-elevation regions, which are already vulnerable to climate change and are currently experiencing its impact, will have severe consequences in the future.
Therefore, biological invasions will add pressure and increase risks to the most vulnerable ecosystems in Nepal.
Along with climate change, anthropogenic disturbance is considered a major driver that promotes plant invasion in mountain ecosystems (Davis, Grime, & Thompson, 2014;Pauchard et al., 2016). Nepal has undergone a significant transformation due to infrastructure development, tourism and trade expansion (Lennartz, 2018;Nepal, 2000). In the mountain and lowland regions of Nepal, newly built roads destabilize slopes and trigger landslides, creating bare ground suitable for colonization by the IAPs (Lembrechts et al., 2012;Lennartz, 2018). Roads play an important role in the spread of alien species by facilitating dispersal pathways and by providing disturbed sites for percolation from roadsides into the natural adjacent vegetation (McDougall et al., 2018). With tourism industry predicted to be grown in the future, human mobility, trade and transport will increase significantly. This may promote the dispersal of IAPs from lowlands to high-elevation regions and to new areas in the lowlands. All of these socio-economic transformations favour the spread of IAPs, and climate change will open up suitable regions by reducing climatic barriers for them to invade higher elevation zones (Hellmann et al., 2008;Pauchard et al., 2009  . Future climate will increase the distribution of the IAPs that were ranked by local communities as the worst, such as Ageratum houstonianum, Chromolaena odorata, Ageratina adenophora and Mikania micrantha . Furthermore, the distribution of Parthenium hysterophorus and Lantana camera, which are considered the most troublesome weeds in the region , will also increase in the future. Therefore, the economic loss and negative impacts caused by IAPs on food security, livelihood, biodiversity and ecosystem services in the future may be augmented if preventive and control measures are not immediately taken seriously. The models developed from the current distribution are extrapolated in time and space to forecast potential IAP invasions under future climate and may not capture the issue of non-analogous climatic space (Fitzpatrick & Hargrove, 2007). Shifts in species range involve multiple ecological processes such as dispersal, demography, physiology, species interactions, population interactions and evolution operating at multiple scales (Urban et al., 2016). Furthermore, the correlation structure of future climatic conditions could be different from current conditions, thereby leading to errors in predictions.
Therefore, SDMs do not explicitly consider these uncertainties caused by non-analogous climate space and ecological processes that affect the species (Elith & Leathwick, 2010). Limiting the areas within the current extent of occurrence (e.g., a MCH) in analysing the change in climatically suitable niches under current and future climate prevents severe changes in the total amount of suitable area (Wright et al., 2015). Furthermore, there are other potential issues such as modelling algorithm (Elith et al., 1999), the choices of environmental variables used (Synes & Osborne, 2011), for future climate, GCMs used (Steen, Sofaer, Skagen, Ray, & Noon, 2017), collinearity (Dormann et al., 2000), model complexity (Wright et al., 2015), model evaluation method (Lobo et al., 2013) and threshold values to produce binary maps (Liu et al., 2016) that can influence model outcomes. In addition, future land use change scenarios can also alter future species distributions (Martin, Dyck, Dendoncker, & Titeux, 2017). There is no agreement on optimal ecological modelling strategy, and such a strategy is unlikely to emerge due to the context-specific nature of the modelling process (Heikkinen et al., 1982;Wright et al., 2015). Despite the uncertainties, it was argued that some amount of model extrapolation for ecological management in a changing climate is essential for practice (Mahony, Cannon, Wang, & Aitken, 2008). Improvements of models are a crucial issue for enhancing the predictive accuracy of the models.
Despite global and local efforts to manage biological invasions, the number of alien species has been ever increasing across all taxonomic groups and geographic regions of the world (Seebens et al., 2017). Climate change has a potential to create more favourable regions in the future for IAPs as shown by this research and other studies (O'donnell et al., 2011;. By creating climatically suitable regions in the most vulnerable natural and agro-ecosystems that provide essential ecosystem services, climate change is likely to amplify the impacts on ecosystems and economy in the future by two major ways. First, climate change negatively affects ecosystems and native species by changing their distribution, composition and phenology (Walther et al., 2002) and hence reduces their resilience to biological invasions. Second, climate change facilitates the encroachment of invasive species by removing current climatic barriers (Hellmann et al., 2008). Cold temperatures limit invasion by many alien species in high-elevation regions , but climate change will elevate this barrier to a higher elevation. The increase in biological invasions will have a serious consequence on the country's economy and local livelihoods. Given that little attention is paid to biological invasions in biodiversity conservation and climate adaptation policies in Nepal, where the formulation of a national strategy for the management of invasive species is still underway, the result of this study as a precautionary note might be helpful to formulate such policies. With this result, we urge that early detection and preventive actions should focus on the mountainous areas of the country. Apart from distribution modelling, a better understanding of species traits, dispersal pathways and the mechanism of the natural filters that prevent colonization of invasive species, as well as the community perceptions and involvement in management, are necessary. Our results show a diverse response of IAPs to climate change; therefore, species-specific prioritization exercises may be helpful to better manage and monitor specific IAPs.

ACK N OWLED G EM ENTS
We used species occurrence data gathered in research projects awarded to BBS by International Foundation for Science (Sweden, Grant no. 5306-1), National Trust for Nature Conservation (Nepal) and International Centre for Integrated Mountain Development (ICIMOD). In addition, we also used some occurrence data gener- Srees and S Pujara for their assistance in the field to collect data.
We are also thankful to Prithvi Narayan Shrestha (The Open University, UK) and Lily Jane Clarke (University of Montana, Montana, USA) for English language edits. We thank two anonymous reviewers for their useful comments that have helped to improve the manuscript significantly.

DATA ACCE SS I B I LIT Y
The occurrence data of the IAPs of Nepal used for analysis in this article are available at Dryad repository; https ://doi.org/10.5061/ dryad.27257r3.