Modelling Hotspots for Invasive Alien Plants in India

Identification of invasion hotspots that support multiple invasive alien species (IAS) is a pre-requisite for control and management of invasion. However, till recently it remained a methodological challenge to precisely determine such invasive hotspots. We identified the hotspots of alien species invasion in India through Ecological Niche Modelling (ENM) using species occurrence data from the Global Biodiversity Information Facility (GBIF). The predicted area of invasion for selected species were classified into 4 categories based on number of model agreements for a region i.e. high, medium, low and very low. About 49% of the total geographical area of India was predicted to be prone to invasion at moderate to high levels of climatic suitability. The intersection of anthropogenic biomes and ecoregions with the regions of 'high' climatic suitability was classified as hotspot of alien plant invasion. Nineteen of 47 ecoregions of India, harboured such hotspots. Most ecologically sensitive regions of India, including the 'biodiversity hotspots' and coastal regions coincide with invasion hotspots, indicating their vulnerability to alien plant invasion. Besides demonstrating the usefulness of ENM and open source data for IAS management, the present study provides a knowledge base for guiding the formulation of an effective policy and management strategy for controlling the invasive alien species.


Introduction
Invasive alien species (IAS) pose serious environmental problems by disassembling natural communities, and adversely affecting ecosystem structure and function [1][2][3][4]. The negative impacts of IAS on native biodiversity, economy, and human health are well-established [5,6]. Major drivers of global change such as climate warming, deforestation, habitat fragmentation, changes in land use and land cover, rapid economic development, and population explosion accelerate the invasion process [4]. Several different approaches have been taken to the study of alien species in India during the last two decades. These include floristic survey and documentation, ethnobotanical description, and experiments in the disciplines of ecology, ecophysiology and genetics [4,[7][8][9][10][11][12][13][14][15][16][17][18]. However, most of these studies were conducted at a local scale, and are species-specific, providing limited scope to construct a pan-Indian IAS scenario. In general, research efforts on alien species invasion in the Indian sub-continent are far behind in scope, intensity occur. These 295 species were checked for availability of georeferenced global occurrence records for model training, and accordingly 155 species were selected. Species occurrence data from Global Biodiversity Information Facility (GBIF) were used for model building. Of the total 825 global terrestrial ecoregions [38], 47 ecoregions located within the political boundary of India were selected. Similarly, of the 35 global biodiversity hotspots [39], 4 hotspots in India were selected viz., Himalayas, Western Ghats, and parts of Indo-Burma and Sundaland. All the 21 anthropogenic biomes under six groups [40] were included in the analysis for delineation of invasion hotspots.

Geographical extent of analysis
We selected Africa, Australia, Europe, North America and South America as training regions for the niche models (Fig 1). Data from species' entire distribution ranges provide a closer approximation of the fundamental niche [41]. Therefore, occurrence data points spread over a large geographic region e.g., at continental level, represent adequate sampling extent and encompass broad range of environmental conditions in which the species occurs. These two conditions are essential for robust distribution modeling of the IAS [42]. The modelled bioclimatic niches were then projected to India with a total geographical area of 3,166,414 km 2 .

Delineating invasion hotspots
We delineated the hotspots of invasion in the framework of modelled climatic suitability, ecoregional setting, and anthropogenic influences. Identification of climatically suitable areas for each IAS is fundamental to hotspot delineation [22]. An ecoregion comprises a geographically distinct set of natural communities, sharing a majority of the species, ecological dynamics, environmental conditions, and critical ecological interactions important for their long-term persistence [43]. Since these communities function together as a single conservation unit at regional scale and the impact of IAS invasion on them would be very similar, inclusion of ecoregion as a component for hotspot delineation would facilitate conservation planning. Anthropogenic biomes were considered as the third criterion for hotspot delineation because human influences facilitate the invasion process [2]. Anthropogenic biomes are the result of sustained and direct human interaction with ecosystems, and were identified through empirical analysis of global population, land use, and land cover [40]. Thus, the consensus zone of the ecoregions, anthropogenic biomes, and areas with high climatic suitability for multiple alien invasive plant species represented a hotspot.

Bioclimatic data and grain size
Climate plays a major role in determining the distribution of plants [44][45][46]. We used 19 bioclimatic variables relevant to persistence, phenology and physiology of the selected species to model the climatic niche ( Table 1). The bioclimatic variables are derivatives of monthly temperature and precipitation values. These variables with a spatial resolution of 2.5 arc minutes were downloaded from the worldclim website (www.worldclim.org).
Selection of an appropriate grain size for environmental variables and species occurrence data is important for enhancing the prediction accuracy in ENM [47][48][49][50]. We resampled the predictor layers to a spatial resolution of 0.04°( 4km), and ran the models at a relatively coarser scale because: (i) the geographic scale of species distribution was quite large, (ii) occurrence data for a large number of species were used, (iii) number of occurrence data points varied widely among the species, (iv) interpolated climatic variables were used as predictors, (v) there was a possibility of mismatch between species point data and the corresponding environmental raster data, (vi) possibility of spatial autocorrelation in the environmental data, and (viii) constraints of computation time and resources.

Species occurrence data
About two million global distributional records for the selected 155 species were downloaded from the web resource of the Global Biodiversity Information Facility (www.gbif.org) based on: (i) invasiveness in the subcontinent, and (ii) availability of adequate and well-spread geographical coordinates. These records were checked for spatial errors using Diva-GIS [51]. For each species, only one occurrence per 0.04°grid cell was retained and the multiple occurrences per cell were removed using ENM Tools 1.3 [52]. Finally, 239,236 occurrence records were used in model calibration, evaluation and projection (S1 Table). Presence of these species in India was confirmed through sample field surveys, available literature [13,16,17], and from online resources viz., Global Invasive Species Database (www.issg.org) and Global Species website (www.globalspecies.org).

Model calibration and thresholding
The climatic niche of each selected species was modelled using maximum entropy ecological niche modelling software MaxEnt 3.3.3k [42]. MaxEnt estimates the maximum entropy probability distribution function based on the input environmental variables in order to predict the area of potential distribution of a species. It defines the ecological niche boundary of the species by constraining the probability distribution based on the environmental parameters of the grid-cell presence record [42]. The attributes of MaxEnt making it the most widely used software for ecological niche modelling are: a strict mathematical definition, continuous probabilistic output, simultaneous handling of continuous and categorical environmental data, functionality to investigate variable importance through jackknife procedure, capacity to handle low sample sizes, facility to allow cross-validation, bootstrapping and repeated subsampling in order to test model robustness, multivariate environmental suitability surfaces tool to identify novel climate conditions, and simplicity of model interpretation [42,[53][54][55]. The Maxent model was run with default settings for the following parameters: number of background points used (10,000), number of iterations (500), convergence threshold (0.00001), and regularization multiplier (1). We selected the 'auto features' option in Maxent software as it helps in selecting the appropriate model fitting functions in respect of each of the species. To derive the optimized model, we executed 10 replicated runs for each species by employing a cross validation technique. In cross validation, the occurrence data is randomly split into a number of folds. Models are generated leaving out each fold in turn, and the left out folds are used for intrinsic evaluation [56]. Cross validation technique was followed to optimize the model by executing 10 replicated runs that generated average, maximum, minimum, median and standard deviation. Model quality was evaluated based on AUC value, and the models were graded as poor (AUC<0.8), fair (0.8<AUC<0.9), good (0.9<AUC<0.95), and very good (0.95<AUC<1.0) [57]. The probabilistic models were truncated by applying the 10 percentile training presence logistic threshold rule in MaxEnt, which generated species-specific binary maps of climate suitability and unsuitability [58]. Novel climatic conditions in India compared to the model training continents were identified and analysed using ExDet Tool Ver. 1.0 [54]. ExDet (Extrapolation Detection) is a robust tool which helps in detecting and quantifying the environmental novelties i.e. novel univariate range (Type 1 novelty) and novel combinations of covariates (Type 2 novelty), in the projected landscapes.

Multispecies climatic suitability maps and identification of invasion hotspots
Multispecies climatic suitability maps were generated by summing up fair, good, and very good quality thresholded models in ArcGIS platform. These climate suitability maps comprised six different sets of consensus models characterized by different degrees of model agreement. In a given consensus model, areas having maximum model agreements were considered suitable for occupancy by maximum number of IAS. In total six consensus maps were generated, of which five maps depicted the projections for Africa, Australia, Europe, North America and South America, and the sixth map depicted a generalized picture of climate suitability for IAS. The sixth map was generated by summing up the multispecies projections of all five model training continents. The species-specific models for each of the 5 continents (as per availability), that represented a 'climate information bit' were combined to obtain the climatic niche of all 155 IAS. Summing-up of these information bits in the projected landscape yielded a consensus map of 525 models with a range of model agreements (15-284 model combinations) which were subsequently categorized as high, medium, low and very low through quantile classification in ArcGIS. The region with the maximum agreements was designated as 'high' category, and was deemed to be suitable for a large number of IAS. The alternate technique for the above classification could have been to use global occurrence data of each species for species-wise projection models and subsequent generation of a consensus map. However, given the fact that most of the species selected were not uniformly distributed throughout the world, such an approach would have resulted in environmental bias. This bias arises out of the spatial bias caused due to collection of non-random occurrence data and background sampling [59].
The classified climatic suitability maps were then overlaid on the terrestrial ecoregions [43], (The Nature Conservancy, http://maps.tnc.org/gis data.html), and the anthropogenic biomes [40], (CIESIN: http://sedac.ciesin.columbia.edu). The ecoregions where >50 percent of the area is climatically suitable with a 'high' consensus level, and contain at least three anthropogenic biome types, were considered as 'invasion hotspots'. The identified invasion hotspots were overlaid on biodiversity hotspot shape files (available at www.conservation.org) using ArcGIS software, to identify the areas vulnerable to alien species invasion. A KML version of the invasion hotspots was created using DIVA GIS 7.5 [51], which was then superimposed on Google Earth 7.1.2 to identify vulnerable forest areas.

Model performance
Based on the AUC values, species-specific models from each of the continents were graded as very good, good, fair, and poor. More than 85 percent of the AUC train and AUC test values for all the continents occurred in the AUC range of 0.8 to 1.0 that discriminated the suitable and unsuitable climates at fair, good and very good levels (Fig 2, Table 2).

Climate novelty
The ExDet tool detected maximum climatic novelty (Type 1) in India for bioclimatic background projections of Europe and South America, while relatively less novelty was detected in the projections from Africa, Australia and North America (S1 Fig). Though most regions in India have 'Type 1' climatic novelty, novel combination of bioclimatic variables i.e. Type 2 novelty, was detected in some areas of Western Ghats, northeastern India and Western Himalaya (S1 Fig). The climatic novelty in these area was due to temperature and precipitation seasonality (S2 Fig).

Climate suitability and invasion hotspots
Continent-specific niche model projections depicted varied pictures of climatic suitability in India. Projections from Africa and Australia mostly predicted the northeastern region and eastern coasts of Peninsular India to have 'high' climate suitability. The North and South American projections matched with the peninsular and northeastern region of India. Western Himalaya was predicted to have 'high' climate suitability with the projections from South America. However, projections from Europe presented a slightly discordant picture with areas of 'high' climatic suitability in the Western Himalaya, Terai areas bordering India and Nepal, Aravalli range in Rajasthan, Hills in the Eastern Ghats bordering Andhra Pradesh and Odisha, and the Naga and Khasi Hills in northeastern India (Fig 3).
The combined projection from all the 5 continents revealed that the coastal areas, the northeastern region, and the Western Himalaya have 'high' climatic suitability for diverse IAS. The projections reveal that more than half of Andaman and Nicobar Islands, Andhra Pradesh, Assam, Dadra and Nagar Haveli, Daman and Diu, Goa, Kerala, Manipur, Meghalaya Mizoram, Nagaland, Odisha, Pondicherry, Tamil Nadu, Tripura and West Bengal, have 'high' climatic suitability for IAS (Fig 4).
Overlaying the combined projection on the ecoregions and major biome types in India revealed their potential to harbor diverse IAS (S2 Table,  The identified invasion hotspots coincide the major biodiversity hotspots of India i.e. Western Ghats, Indo-Burma, and Eastern Himalaya (Fig 5). They contain diverse signature of human activities, and are characterized by 5 major anthropogenic biome types and 18 subtypes (Fig 6). Apart from the biodiversity hotspots, the invasion hotspots also coincide forest reserves, islands, coastal forests and mangrove ecosystems of India.

Discussion
By using a climate-based approach, recent research has identified species having invasion potential in climatically similar regions [20,22,24,31]. A synthesis of data on the characteristics of invasive species and identification of correlates of invasion within and across different biological groups also revealed that climate/habitat match is the only characteristic that is significantly and consistently associated with invasive behaviour across biological groups including plants [37]. While the utility of climatic niche at large spatial scales such as regional, continental and global is well-established, the role of other non-climatic factors in determining the presence or absence of a species at finer i.e. regional or local scales is also supported by many studies [33,34]. Climate principally affects the distribution of invasive species as they can only succeed in regions where they are not limited by climatic constraints. Moreover, the conservatism of climatic niche makes the climate-based model predictions widely applicable [53]. Considering the above facts and in view of the major advantage of ENM in assessing the invasion potential of large numbers of species, in the present study, climatic niches of the multiple IAS were modelled [60].
Comprehensive distribution maps depicting climatic suitability and invasion hotspots in India were generated through niche models of a broad spectrum of species covering different taxa, life-forms, and plant functional types. The continent-specific multispecies ENM projections show some discordance in the scenarios of climatic suitability for IAS in India (Fig 3). Such discordance may principally be attributed to variation in the ecological niche models resulting from the inbuilt correlation structure of the bioclimatic datasets of each continent [61]. Acknowledging the different sources of origin [13,15,16], building such scenarios would help us in having a better understanding of alien species invasion. Earlier workers reported that 35 percent of the alien plants in India are from South America, and the rest are from Asia (21%), Africa (20%), Europe (11%), Australia (8%) and North America (4%) [17].  Nevertheless, the continental projections showed some consensus areas that helped us in identifying the invasion hotspots in India. The model consensus clearly outlines the following:

Invasion hotspots coincide with biodiversity hotspots
The concurrence of invasion hotspots and biodiversity hotspots may be attributed to similar environmental requirements of the alien species to that of their native counterparts [62].  Nevertheless, presence of an environmental gradient providing variability in resource availability may promote species invasion in the biodiversity hotspots [63][64][65].
The biodiversity hotspot regions in India have the highest concentration of the country's forest cover and include most of the protected areas [66]. In the last three decades, these biodiversity hotspots have experienced sizeable forest cover loss [67] because of biotic pressure, shifting cultivation, agricultural expansion, and urbanization [68]. Such disturbances create a state of disequilibrium in the environmental conditions, which make these areas vulnerable to invasion [69][70][71][72][73]. Such vulnerability of ecosystems to invasion may be explained based on the concept of availability of 'vacant niches' i.e. unfilled opportunities for additional species in ecosystems or by creating new 'ecological opportunities' for occupancy by the IAS [74,75]. 'Ecological opportunities' denotes a pulse in the availability of readily assimilable nutrients in the novel habitats, essential for colonization, growth and development of IAS [76]. In fact, concurrence of invasion hotspots and the biodiversity hotspots is an indication of impending threat posed by the IAS to the global biodiversity hotspots. Human activities have been argued to facilitate species invasion [77,78]. The identified invasion hotspots in this study also coincide areas with intense anthropogenic activities, which is evident from the dominance of anthropogenic biomes such as croplands, rangelands, urban areas, and rainfed villages (Fig 6). Thus, anthropogenic activities seem to be a strong correlate of invasion hotspots.

Invasion hotspots coincide with islands, coastal forests, freshwater swamp forests, mangroves, and forest reserves
The invasion hotspots coincide with the Andaman Islands, the coastal forests, swamp forests, the mangrove ecosystems, and numerous forest reserves in India. The explanation for concurrence of these ecosystems with the invasion hotspots may be the same as mentioned above [62]. However, this concurrence should be viewed as an imminent threat to the native biodiversity of the islands, coastal forests, swamp forests, mangrove ecosystems, and forest reserves. The vulnerability of island ecosystems to alien species invasion has been attributed to the following reasons: (i) Islands have high net resource availability, and the native species have a poor ability to preempt those resources [76]. Consequently, the competitively superior invaders get an opportunity to occupy the vacant niches. (ii) Islands are exposed to high propagule pressure and the native species face a strong inter-specific competition from the alien invasive species [79], which make them susceptible to invasion. By and large, the above ecosystems are under constant threat of alien plant invasion due to landuse changes and biotic factors [80][81][82]. For example, although mangrove ecosystems in general are not susceptible to invasion, the mangrove forest areas adjoining terrestrial ecosystems, at the heads of estuaries, and areas having open canopy cover due to anthropogenic disturbances are susceptible to invasion. Some of the major port towns such as Mumbai, Ratnagiri, Panaji, Nagapattnam, Chennai, Kakinada, Paradip, and Haldia, fall within the identified invasion hotspots. Such areas may provide suitable habitats for colonization of IAS after being introduced through shipping routes [83]. Besides, the invasion hotspots also coincide with many hill and forest ranges, reserve forest areas, and protected areas in the Eastern Ghats, which have high floristic diversity. Some of the important forest areas and protected areas include:

Policy implications
Our work has illustrated the usefulness of open source data and software in generating a knowledge base having practical implications of IAS management by understanding and analyzing the alien species invasion problem which operates at a broad geographical scale. Some of the inferences that have strong policy implications are: 1. Invasion hotspots coincide with hotspots of biodiversity: This can be viewed from a threat perspective as it would endanger the persistence of many resident and endemic species. Enumeration of IAS should be undertaken in the biodiversity rich biomes and their impact on other native species at ecosystem level needs to be analysed.
2. Invasion hotspots coincide with islands, coastal areas and protected forest areas: Surveillance and quarantine measures should be stepped up in such areas as they may act as important sites of introduction, colonization and establishment of IAS.
3. Invasion hotspots coincide with croplands, rangelands and village biomes: Efforts for eradication and control of the alien invasives in these areas should be given top priority as it is relatively easier to control them at the beginning of their spread.
4. The remote forest biomes are unlikely to be invasion hotspots: The remote forest biomes did not coincide with any invasion hotspot. These biomes in India include areas with high tree cover and low human population. Therefore, it is important that such biomes are maintained in its present condition with minimum anthropogenic interferences.

Caveats of the study
While recommending the above, we do recognize the following limitations of our predictions on hotspot of alien species: (i) although a biome may not support multiple IAS, yet a particular invasive species can pose serious threat to the biome, and (ii) the predictions are scale-dependent i.e. applicable to regional and global scales. Extrapolating invasive species distribution models to a different geographical region relies on two fundamental assumptions of ENM: (i) species are at equilibrium with their surrounding environment, and (ii) they tend to retain their native range niche preferences in the invaded range [84]. Species are considered to be at equilibrium with the climate if they occur in the entire breadth of the climatically suitable areas [85]. Although niche models generally assume that the distributions of the species are at equilibrium with current climate [27], several workers have shown that the validity of this assumption varies substantially across different groups of organisms [86]. In the present study, most of the selected species are likely to be closer to equilibrium. Although empirical data for supporting this is not available, likelihood of the IAS attaining equilibrium with the climate can be argued based on the following evidences: 1. Minimum residence time: Minimum residence time of a species is the time since its first arrival in the invaded range. A longer residence period would give sufficient time for colonization and spread, and reach equilibrium with climate. It is usually ascertained through: (a) published literature, and (b) the date of the first herbarium specimen collected for the species [87]. Although data on the exact time of introduction for IAS is scarce, there is a strong correlation between known time of introduction and the date of the first herbarium collection [88]. In the present study, the residence time of 86 species was determined as more than 150 years based on the published literature and herbarium records [89][90][91][92][93][94]. Their introduction to India can be linked with the historical trade relationship with Britain, Portugal, Spain, France, Middle East and Central Asia [17]. Irrespective of the residence time, IAS do vary greatly in the rate of their spread based on species-specific traits. For instance, at least two species out of 155 viz., Alternanthera philoxeroides and Mikania micrantha which were introduced to India only in 1940s, have become highly invasive and have spread throughout the country [4,95].
2. Species traits: Life-form and life-cycle duration e.g., annuals, which reproduce each year will have more generations than trees, even if they have arrived at the same time. At least 79 IAS out of 155 are annual herbs with r-strategic life-history having ability to produce large number of propagules, with short life span and with multiple generations within a limited time period. The extent of occurrence of each of the modelled species has a wide native range, thus having high invasive potential [15,96].
3. Mode and rate of dispersal: Mode and rate of dispersal of the propagules determine how fast the IAS would reach and colonize new areas [97]. Inclusion of proxy variables representing various processes that contribute to species spread viz., network of roads has been shown to enhance model accuracy [98][99][100]. The high density of road network and human habitations in India is likely to facilitate localized dispersal of IAS in the vicinity of occupied sites, besides contributing towards high propagule pressure. The road network and human habitations also enhance the invasion potential by creating disturbance gradients at landscape level [101]. Further, the spatial extent and rate of dispersal of propagules of selected IAS are likely to be augmented by biotic agents, the predominant mode of dispersal, especially the highly mobile human beings [102].
Thus, great dispersal abilities coupled with a reasonably longer residence time of the selected suite of IAS, reflects its high likelihood of being in equilibrium with the climate.
Responses of plant species, in general, to climate change are more likely to be accurately forecasted by models correlating present day distributions with climate than the animal groups [86]. IAS benefit from a variety of global change factors including enhanced nitrogen deposition, increased atmospheric CO 2 concentrations, habitat fragmentation and other disturbances [103][104][105]. In addition, as argued by Colautti and Lau [106], the contemporary evolution may provide invasive species an evolutionary advantage over native species because of: (i) better adaptation due to co-evolution of IAS and humans [107], (ii) the act of multiple introductions from different geographical sources that increase genetic variation in physical and ecological traits, and produce novel phenotypes facilitating adaptation in introduced regions [108], (iii) greater population size and high density of IAS enhances the probability of adaptive mutations against changing environment ensuring competitive advantage over natives [109], and (iv) invaders escape many evolutionary constraints by leaving behind their specialized natural enemies in the native range [110], providing selection opportunity on multiple traits and more freedom to adapt in comparison to the native species.
In the present study, the ExDet analysis shows certain areas in the Western Ghats, Western Himalaya, and northeastern region of India as novel in some climatic aspects (S1 and S2 Figs). The predicted climatic suitability for IAS in these regions, and their presence as validated through field surveys indicates the robustness of the predicted hotspots. However, to have further explicit and mechanistic explanations, a fine grain empirical analysis at ecosystem, provenance and intraspecific level is required [111].
The above evidence suggests of possibilities, where species evolve, often rapidly under selection, leading to differences in population traits from original introduction [112,113]. Further, local adaptations also enable the IAS to invade novel regions [111] resulting in the expansion of distributional range. These processes are over and above the likelihood of being in equilibrium with current climatic conditions, and thus would have little bearing on the current study.
The availability of empirical data pertaining to these variables and processes in future would further enhance model accuracy and predictions.
Supporting Information S1 Fig. Extrapolation detection (ExDet) analysis showing areas with two types of novel conditions in the Indian subcontinental region against the model training regions: Africa (A), Australia (B), Europe (C), North America (D), and South America (E). Regions with red color have 'Type 1' novelty, and blue color have 'Type 2' novelty. Regions in green color represents similar environmental conditions. The 'Type 1' novelty represent points that are outside the range of individual covariates, and 'Type 2' novelty represent the points that are within the individual covariate range but have novel combinations between the covariates. (TIF) S2 Fig. Extrapolation detection (ExDet) output map depicting the distribution of most influential covariate (MIC) out of the 19 bioclimatic variables which contributed to the 'Type 1' and 'Type 2' novelty in the Indian subcontinental region against the model training regions of Africa (Fig A), Australia (Fig B), Europe (Fig C), North America (Fig D), and South America (Fig E). Individual colours along with numbers below the colour ramp depict the MIC i.e. the bioclimatic variables. Blue areas in the map, labelled as 'none' in the colour legend, do not have any covariate which falls outside the range of native range bioclimatic data or have any non-analogous combination of these covariates.