Ageratina adenophora and Lantana camara in Kailash Sacred Landscape, India: Current distribution and future climatic scenarios through modeling

The Himalayan region is one of the global biodiversity hotspots. However, its biodiversity and ecosystems are threatened due to abiotic and biotic drivers. One of the major biotic threats to biodiversity in this region is the rapid spread of Invasive Alien Species (IAS). Natural forests and grasslands are increasingly getting infested by IAS affecting regeneration of native species and decline in availability of bio-resources. Assessing the current status of IAS and prediction of their future spread would be vital for evolving specific species management interventions. Keeping this in view, we conducted an in-depth study on two IASs, viz., Ageratina adenophora and Lantana camara in the Indian part of Kailash Sacred Landscape (KSL), Western Himalaya. Intensive field surveys were conducted to collect the presence of A. adenophora (n = 567) and L. camara (n = 120) along an altitudinal gradient between 300 and 3000 m a.s.l. We performed Principal Component Analysis to nullify the multi-colinearity effects of the environmental predictors following MaxEnt species distribution model in the current and future climatic scenarios for both the species. All current and future model precision (i.e., Area Under the Curve; AUC) for both species was higher than 0.81. It is predicted that under the current rate of climate change and higher emission (i.e., RCP 8.5 pathway), A. adenophora will spread 45.3% more than its current distribution and is likely to reach up to 3029 m a.s.l., whereas, L. camara will spread 29.8% more than its current distribution range and likely to reach up to 3018 m a.s.l. Our results will help in future conservation planning and participatory management of forests and grasslands in the Kailash Sacred Landscape–India.


Introduction
Invasive alien species (IAS) rank among the top three threats to global biodiversity, the other two being unsustainable harvest of various species from the wild and habitat degradation and

PLOS ONE
Sanctuary) which provides important ecosystem services for the region. However, the non-PA category of forests includes reserve forests and protected forests under the control of the State Forest Department of Uttarakhand, managed for an extended period for the production of timber and non-timber forest products (NTFPs). KSL-India also includes highly interspersed small forest fragments falling under two broad categories, i.e., civil/soyam forests and community forests are under the control of the Revenue Department of the State Government. The latter category of community forests, for nearly nine decades or so, has been managed by the local communities and is referred to as 'Community Forests' or 'Van Panchayat Forests'. This community forest serves as the primary source of livelihood for locals. The KSL-India has experienced a considerable change during recent decades in terms of land use and land cover, along with the development of infrastructure and patterns of natural resources used by the local communities [23,24].

Species presence
After the extensive vegetation survey and consultation with experts, the two most widely spread IAS, i.e., Ageratina adenophora (Asteraceae) and Lantana camara (Verbenaceae) were selected for the present study. A. adenophora is native of Central America [25] and distributed worldwide [26]. L. camara, native to Tropical America [27] is considered to be one of the most troublesome invasive species worldwide [28]. A. Adenophora, a perennial, herb produces a large number of tiny seeds which are dispersed by wind and human activities [29]. Its ubiquitous property to invade diverse habitats, such as roadsides, riverbanks, forest edges, crop fields, wastelands, and rubbish dump edges makes it fit to invade various ecosystems. L. camara is one of the most widely distributed IAS in India [30], resistant to fire and can regrow if burnt [31]. It reproduces primarily by seeds as well as through coppice and prefers to grow in degraded habitats. Though, the species selected for the study are known to negatively affect the wildlife habitats and native vegetation, they are also known to have a few minor uses to local communities. Leaves of A. adenophora are used for cattle bed [32], or its leaf paste is applied to cuts and wounds and leaves are also have the properties of antimicrobial activity, antipyretic activity, wound healing, antioxidant activity, analgesic activity, anti-tumor activity, insecticidal activity, antiviral activity, anti-inflammatory activity, and larvicidal activity [33,34]. The wood of L. camara is used as fuel [35,36], a source of food for birds [37][38][39] and nectar for butterflies and moths [40,41] as flowering and fruiting happens throughout the year.

Presence record
We recorded A. adenophora and L. camara presence between 300 and 2500 m a.s.l. covering a wide range of habitats including agricultural lands, grasslands, fallow lands and forests, proximity to major road networks and naturally occurring drainages in different parts of the landscape during 2014 to 2017. A total of 567 and 120 presence locations for A. adenophora and L. camara, respectively were recorded using handheld Garmin N72 GPS (Supplementary Information). Of these, 80% locations were used for probability distribution modeling and the rest 20% were used to evaluate the model performance. The study was conducted outside the protected area and no Schedule-I or -II species were used in the study.

Selection of environmental and bioclimatic variables
We used two-time periods ("current" and "2050") environmental data to model the A. adenophora and L. camara current and future distribution in KSL-India. A total of 22 environmental variables were used for probability distribution modeling in each period. We retrieved standard 19 bioclimatic variables for WorldClim version 2 [42]. These are the average for the years 1970-2000 with a spatial resolution of 30 seconds (~1 km 2 ). These variables are considered as current bioclimatic variables. We also retrieved 19 bioclimatic variables with a spatial resolution of 30 seconds (~1 km 2 ) for the year 2050 (average for 2041-2060). These are the climate projections as per the Intergovernmental Panel on Climate Change (IPCC) 5 th assessment report from Global Climate Models (GCMs) for one Representative Concentration Pathways (RCPs). However, we used RCP 8.5. Because, out of several IPCC climate change scenarios, RCP 8.5, generally taken as the basis for worst-case climate change scenario among all (viz., RCP 2.6, RCP 4.5, RCP 6.0 and RCP 8.5.) [43][44][45]. The data were downscaled and calibrated using WorldClim 2 as baseline 'current' climate [42,46]. These bioclimatic data are also the most recent GCM climate projections that have already been considered in the IPCC fifth assessment report.
We also considered four major static topographic factors for both the "current" and "2050" time frame. These are elevation (ASTER, Digital elevation data from Advanced Space-borne Thermal Emission and Reflection Radiometer), slope, aspect, Euclidian distance from water channels and Euclidian distance from major roads. These variables were assumed to be constant across the entire study period. The details of all the variables used in this study are provided in Appendix 1 in S1 Appendix.

Pair-wise Pearson correlations and principal component analysis
We checked for multi co-linearity with pair-wise correlations for both "current" and "2050" environmental data sets to avoid spurious model calibrations [47]. It is common practice to retain only predictors with pair-wise correlation < |0.7| [48]. To nullify the correlation structure of the variables and inter-relationship among them, we implemented the Principal Component Analysis [PCA; 49,50]. The PCA reduces the number of orthogonal predictors each of which is the linear combination of 24 original environmental variables. All the exploratory variables were subjected to PCA to extract the factors of significant contribution using prcomp function and predicting the function of R version 3.3.3 [51]. Furthermore, we obtained 24 Principal Components (PCs) in raster format. For both "current" and "2050" data sets, the first six PCs accounted for 99% of the variability and were chosen for further multivariate species distribution modeling (Table 1).

MaxEnt species distribution modeling
Maximum Entropy (MaxEnt) distribution modeling was implied to identify the current and future invasion for both, A. adenophora and L. camara in KSL-India. MaxEnt is a machine learning algorithm that is widely used for species distribution modeling that uses presenceonly data [52,53] to predict the potential distribution of a particular species in the current as well as future climatic scenarios [54,55].
The occurrence data from the field was collected for both A. adenophora (n = 567) and L. camara (n = 120). One thousand background points were drawn randomly across the KSL-India by considering the pseudo-absence of both species. We generated four main MaxEnt models (two for each species) using the six PCA derived explanatory variables for both the 'current' and '2050' year.
We implemented k-fold cross-validation [56] to build each main MaxEnt model. Firstly, we used standard k-fold cross-validation in our randomly partitioned approach. In our k-fold cross-validation approach, the occurrence localities were divided randomly into 10 bins, where each bin constitutes an equal proportion of sample size [57,58]. Then, each model was developed iteratively, using (k − 1) bins for model calibration in each iteration, with the remaining fold retained for evaluation. This is repeated until all the folds were utilized at least once for model evaluation. MaxEnt produces a model based on a series of 'features', such as linear, hinge, quadratic, product, threshold and discrete; we used all except for discrete, as all of our explanatory variables were continuous [59]. In the 'feature' functions, we set the beta multiplier as 0.5 (medium) that affects the smoothness of the model output. We have opted for the clamping function in each model as the use of the function is strongly suggested for projecting future species distributions [60,61].
As we set the models via k-fold cross-validation (k = 10), for each dataset, we used the respective evaluation localities to calculate Area Under the Curve (AUC) of the Receiver Operation Characteristic (ROC), which is a measure of the overall discriminatory ability of the model or to precisely evaluate model performance by plotting its model sensitivity graph [62]. We averaged those values from each sub-models of each main model and represented the averages for each of the four main MaxEnt models. The AUC score closer to 1 indicates that the model had accurately predicted the habitat, while a value �0.5 indicates the low accuracy of the model [63]. We also calculate other model evaluation matrices such as kappa max (kappa), true positive rate (TPR), true negative rate (TNR), and True skill statistics (TSS) using the 'dismo' package in R [64].
Binary habitat invasion maps for the 'current' and future '2050' were prepared by using the threshold of maximum training sensitivity plus the specificity logistic threshold. The threshold approach of 'Maximum training sensitivity + specificity' is a more reliable, restrictive and conservative approach to understand the potential distribution of a species. This particular threshold approach is one of the popular methods either for presence/absence or presence-only data [65]. We used 'dismo' package in R [64] along with a source file developed to build all MaxEnt models [66] as well as to evaluate their performance.

Invasion across the elevation gradient
To check the future spread of A. adenophora and L. camara, we extracted all elevational range values from a 30 m resolution digital elevation model (ASTER) using the threshold polygons generated from MaxEnt main models. Elevation gradients that fall within the distribution ranges were re-classified into multiple classes considering 100 m a.s.l. intervals. Area invaded by both species in each elevation class was compared between current and 2050 year using paired sample t-test. The distributions of elevation class vs. invaded areas of both A. adenophora and L. camara were plotted separately for the current year and 2050 to visually assess their spread pattern across the elevation classes.

Principal components
The

Results
The pair-wise Pearson correlations were generally higher (r >|0.70|) for most of the variables (Fig 2). We were not able to use these variables directly in MaxEnt as they showed high correlation, and this avoids the overprediction of suitability in MaxEnt models. However, PCA transformed the 24 environmental variables (each for 'current' and '2050') into a composite set of 24 components that are orthogonal and un-correlated to each other. Six out of 24 components demonstrated the cumulative proportion of variance > 0.99 for both the 'current' and the year '2050' ( Table 1). The first six PCs each from the 'current' and '2050' were represented in Appendices 2-3 in S1 Appendix. Both training and test area under ROC curve (AUC) values were > 0.81 for both species under the 'current' and future '2050' climatic conditions (Fig 3, Table 2). Other model evaluation statistics were provided in Appendix 4 in S1 Appendix. The current occurrence locations of A. adenophora and L. camara were 95.1% and 92.5%, respectively falls within their most climatically suitable current distribution ranges. The relative percent contribution of all PCs from both the 'current' and the year '2050' are presented in Table 3, and all response curves associated with the above MaxEnt model predictions were given in Appendix 5 in S1 Appendix. The current and future distributions of the A. adenophora and L. camara along the elevation gradient of KSL-India are determined by the RCP 8.5 pathway projections for 2050 ( Fig  4A and 4B).
The current distribution of A. adenophora ranges between 212-2616 m a.s.l., while L. camara distributed mostly below 2045 m a.s.l. The potential area that is currently vulnerable to invasion by A. adenophora and L. camara extends to 1403.7 km 2 and 1053.1 km 2 , respectively. If the climate change continues to follow the higher emission RCP 8.5 pathway, by 2050, invasion of A. adenophora would reach up to 3029 m a.s.l., whereas, L. camara will invade up to 3018 m a.s.l. Both A. adenophora and L. camara will spread across the landscape occupying an area of 2039.7 km 2 and 1366.5 km 2 , respectively. The area invasion by A. adenophora (t =

Discussion
Plant invasion is the emerging area of science that impacts substantial economic and ecological imbalance by altering native plant community composition [67] and depleting native plant species diversity [68] thus affecting ecosystem processes [69][70][71]. The lower and mid ecosystems in the mountains are more vulnerable to invasion by exotic plant species due to the altering climate change and increasing anthropogenic pressure [72]. These plant species also exceeded their ranges and many of them have become abundant at higher elevations [73,74]. The increasing incidence of invasion in the high-altitude ecosystems [9] poses a significant threat to the native biological diversity of the region and it is expected to expand in the future [21]. Using principal components as input variables in MaxEnt enhances the prediction of good habitat suitability [50,75]. Alien species that seek to establish themselves in mountains easily proliferate and become difficult to sustain and manage due to the rugged and steppe mountain terrains [76]. It is important to detect potential invaders and take preventative action, since this is the most cost-effective way to restrict the spread of IAPS to new locations. This demonstrated an increasing (both vertical and horizontal) trend of plant invasion in KSL-India. Altering climatic condition and mostly increased temperature has been providing favourable conditions for the invasion and spread of these species in the KSL-India.  The presence of L. camara was recorded widely in forest and scrublands [77,78]. Looking forward, L. camara might also affect the regeneration of native biodiversity in the area. It appeared that in the future both of these IAS would colonize more vigorously in the landscape where currently they are not present. This invasion of both species is likely to happen if climate change continues to follow the higher emission RCP 8.5 pathway. Furthermore, L. camara will spread within a similar altitudinal range where it is currently absent. These species have already spread in mountains and along rivers and roadside [79]. A. adenophora presence between 2800 to 3000 m a.s.l. [13], which suggest that the spread of this species in higher regions of western Himalaya can cause a potential threat to the native biodiversity in the alpine and subalpine zone.
Our study also showed a similar trend, as this species will reach 968 m a.s.l. higher than its current distribution range. Our predictive models also depicted high infestation and expansion in the present climatic scenario by A. adenophora and L. camara in the lower parts of KSL-  India in western Himalaya. The MaxEnt model also predicted that the landscape provides more suitable conditions for the infestation of A. adenophora as compared to L. camara even at its current climatic condition. Both the species will expand their habitats with the increase in temperature and human activity in the landscape. Our study demonstrated similar results as compared to previous studies [80][81][82], but in a more robust and spatially explicit manner. Our study has come up with a current pattern of distribution and established future scenarios of invasion by two highly obnoxious IAS within KSL-India. The distribution maps generated in this exercise would be of much help in further conservation planning and participatory management of forest/grassland management. The findings of our study can be used in the identification of vulnerable areas and important localities which are likely to be infested by these species. Habitat restoration and community-based monitoring at a local scale involving Biodiversity Management Committee (BMC) would be one of the logical steps forward.