Predicting the Distribution of Parthenium Weed ( Parthenium Hysterophorus L .) Under Current and Future Climatic Conditions in Bhutan

Invasion, spread and establishment of invasive alien species in a new environment causes serious ecosystem perturbation and native species extinctions. The problem is further aggravated under climate change as some invasive alien species perform better under elevated temperature and carbon dioxide regimes. Currently unsuitable regions such as high-altitude areas and mountains are likely to become more suited to invasion under future climatic conditions. We have modelled the distribution of one of the most invasive alien species, parthenium weed ( Parthenium hysterophorus L .) which is rapidly colonizing different parts of Bhutan. The study was implemented in R environment using BIOMOD2 package; an ensemble modelling platform for species distribution modelling. Under current climate scenario, about 2.83% (1,099.01 km 2 ) of the country’s total area was predicted to be suitable for parthenium weed invasion, covering 17 out of the 20 districts in Bhutan. Under future climate scenarios, the highest suitability was predicted under RCP4.5 2050 period with about 5,419.69 km 2 anticipated to be suitable. Except for Bumthang, all districts showed suitability to invasions under future climate scenarios. Generally, districts located in the west and south showed more suitability than those in the east and central region. The highest elevation of suitability was predicted to be at 2,931 m above sea level; an upward shift of about 753 m. Based on these findings, there is an urgent need to develop management programs and raise public awareness on the adverse impacts of parthenium weed in Bhutan.


Introduction
The introduction of alien species outside of their native range deliberately or accidentally destabilizes ecosystem functioning [1]. Native species extinctions have been well-documented [2] and debate on whether alien species invasion could lead to mass extinctions of native species is growing [3]. Even the Intergovernmental Platform for Biodiversity and Ecosystem Services (IPBES) of the United Nations has recognized the threats posed by biotic invaders [4]. Further, climate change is escalating the situation by removing climate barriers [5], facilitating the spread and establishment [6] and shifting ranges [7]. The second most threat after habitat loss to biodiversity is Invasive Alien Species (IAS) [8]. Bhutan, a tiny Himalayan democratic constitutional monarchy with an area of 38,394 km 2 , a population of about 0.72 million [9] endowed with very rich biodiversity is no exception and considered spared from invasion by IAS. There is already a number of IAS reported in Bhutan [10][11][12][13][14][15][16] Between 1974 to 2000, from a total of 127 species imported into the country, four species of trees and shrubs and 18 species of grasses have become very serious weeds [15]. There are already about 43 Invasive Alien Plant Species (IAPS) recorded in Bhutan [17] and 16 are a major IAPS [16]. The most comprehensive inventory of IAPS undertaken by [10] reported about 101 species of IAPS in the country. The detrimental effects of IAPS to public health, agriculture, fisheries, and the economy is well-known with losses estimated in millions of dollars. For instance, IAPS such as parthenium weed is known to cause allergies, asthma and bronchitis, and in some cases even death [18]. Consumption of Crofton weed has led to the death of 30 horses [19] and farmers in Bhutan have reported impacts of IAPS on highland pastures [20] and on crops and forests [21]. In Australia, it is estimated that invasive weeds alone cause crop production loss to the tune of AU$ 1.271 billion per year [22], while in Nepal it is estimated that IAPS cause losses about US$ 1.4 billion per year [23] and India, about US$ 91 billion per year [24]. One of the most invasive alien species known to the world with its origin in central Mexico [25] and reported in 45 countries, including Bhutan [26,27] is the parthenium weed. In Bhutan, it was first reported in 1992 along roadsides, in crop fields and settlement areas [28] and it is now reported at elevation as high as 2,600 m Above Sea Level (ASL) [13]. To locals in eastern Bhutan, it is known as "Amartala ngon" (Weed from the place in India called Amartala/Agartala) or "Gungthri ngon" (seeds fallen from the sky) (S. Wangdi, personal communication, 12 March 2021). It reduces crop yield by as much as 50% [27] and causes detrimental effects to public health, agriculture, forestry, and livestock. When used as a livestock feed, milk becomes bitter [29] and meat can become tainted. Under climate change, it is likely to spread because of the benefit accrued from global warming and atmospheric CO2 enrichment [30]. Many invasive species can expand their range under climate change [31,23]. Therefore, there is a need to accurately predict the potential distribution and understand the impact of parthenium weed in Bhutan. This will not only help to make informed decisions on matter concerning biodiversity, public health, agriculture and the economy but also help in early detection to enable prompt actions in order to reduce management costs, which otherwise would become expensive or impossible at later stages. Other than the study conducted by [32] and [33], no study had been conducted on IAPS distribution in Bhutan. Both the studies implemented the most commonly used species distribution modelling program (-MaxEnt) [34] and provided a general overview of parthenium weed distribution at the national level and West-Central region, respectively under current and a future climate (predicted for 2070). To develop management programs and inform policy-decision, further details such as areas under high risk of invasion or currently invaded, percent invasion under current climate and future climates, elevational shifts and distribution predicted for short-(2030s) and medium-term (2050) periods is necessary. Using the BIO-MOD2 package [35,36] and ensemble approach in R v 4.0.3 [37], we have modelled the distribution of parthenium weed in Bhutan for the current (1979-2013) and two future climate periods (2050 and 2070). Additionally, we have analyzed the distribution of parthenium weed at the district level, which will be important while planning for early detection, prevention and control and management programs.

Spatial thinning of occurrence points
The occurrence points of parthenium weed were collected from the three national highways: Thimphu to Sarpang national highway, Thimphu to Punakha national highway and Thimphu to Samtse national highway via the Phuentsholing to Samtse internal road ( Figure 1). The occurrence points were collected using a handheld GPS with 3 m accuracy. The total distance surveyed was about 360 km. The records were also taken for the dried and dead remains of the parthenium weed from the previous years in case the data collection time coincided with the early emergence time at that site. A total of 151 occurrence points were collected between April to August 2020. As spatial autocorrelation can lead to model over fitting and can reduce model performance, we tested the data for spatial autocorrelation and thinned in ArcGIS 10.4 using SDMTOOLBOX v2.4 [38]. A total of 55 occurrence points, which were retained after spatial thinning from 151 records were finally used for modelling.

Selection of predictor variables and testing for multicollinearity
The likely distribution of a species is commonly predicted by using bioclimatic variables and the WORLDCLIM 2 (http://worldclim.org) [39] is the source mostly used to obtain these variables. We have used the CHELSA (Climatologies at High Resolution for the Earth's Land Surface Areas) website (https://chelsa-climate.org/) [40] to obtain the 19 bioclimatic variables; only being that it has out-performed the WORLDCLIM 2 in terms of precipitation, which was demonstrated from a preliminary study using Bhutan [40]. The 19 bioclimatic variables at 30 arc sec (ca. 1 km 2 ) resolution were downloaded for the current (1979-2013), and two future climate periods, e.g. for 2050 (2041-2060 average) and 2070 (2061-2080 average) for Representative Concentration Pathway (RCP) 4.5 and RCP8.5. The RCP4.5 assumes a steady increase in radiative forcing with a projected global mean surface temperature of 1.4°C to 1.8°C, while the RCP8.5 assumes the highest increase in radiative forcing with projected global mean surface temperature of 2.0°C to 3.7°C [41]. The choice of RCPs was also in line with the climate change projection for Bhutan  [42]. Given that Bhutan is a mountainous country with high microclimate variability, we have used the global circulation model, MIROC5 has used by others working in the Himalayan region [43][44][45][46][47]. Additionally, we have generated two topographic variables (slope and aspect) from the Digital Elevation Model (DEM), which was obtained from the Shuttle Radar Topographic Mission (SRTM). Predictor variables were tested for multicollinearity [48] using the vifcor function in the usdm package [49] executed in R v 4.0.3 [37]. There is no consensus as to what an acceptable threshold is, and how it could be used for discarding correlated variables. Based on the recommendation by [50] we have used a threshold of r (Pearson's correlation) >0.8 to screen out the correlated variables. The variable retained after multicollinearity test were: The study did not include other anthropogenic variables such as roads, land use and settlements as future projections of these variables were not available.

Modelling approach and calibration
There are many species distribution models used to model species distributions in a geographic space and environment [51]. We have used an ensemble modelling approach, as it accounts for uncertainties of different models and reported to perform better than a single model [52,36]. The package BIOMOD (BIOMOD2 package in R) developed by [35] has been especially designed to address this issue. The package provides 10 models but for our study, we have only used five, viz., General Additive Model (GAM), Generalized Linear Model (GLM), Generalized Boosted Models (GBM), maximum entropy (MAX-ENT.Phillips) and Random Forest (RF). These models are the most used and are known to perform well [53,54]. The above models in BIOMOD2 were run with these settings. While formatting, the number of pseudo-absences generation was repeated three times, with their number generated set to 2000. Keeping a buffer of 10 km from the presence points, pseudo-absences were randomly selected [55]. Pseudo-absences were required as the models used were based on presence-absence data and we did not have absence data. At the modelling step, we calibrated the models by splitting the data 70:30 (70% for model calibration/training and 30% for model evaluation/testing) [56]. In total, we built 75 models (three pseudo-absence runs, five models and five evaluation runs). In order to build ensemble models, models below True Skills Sta-tistics (TSS) value of 0.6 and Area Under the Curve (AUC) 0.8 were screened out and weights were assigned to each model balanced to their evaluation metric score; what's called a weighted-mean approach. Using TSS scores as a cut-off value, we produced binary maps of suitable and unsuitable areas. Change in suitability areas under future climatic conditions was produced by using BIOMOD2 range size function. All the maps produced in R environment using BIOMOD2 were imported into ArcGIS 10.4 for analysis and visualization.

Model evaluation
We have used the AUC of the Receiver Operating Characteristics (ROC) [57] and TSS to evaluate the model performance. The AUC is a measure of the ratio between the observed presence-absence values and the model predictions and the values range from 0 to 1. This is a threshold independent measure and according to the scores, a model is classified as excellent (0.9-1.0), good (0.8-0.9), fair (0.7-0.8), and poor (0.6-0.7) [58,59]. TSS is a threshold dependent measure and its values range from -1 to +1. Based on TSS, a model is categorized as no better than random (< = 0), poor (0.2-0.5), useful (0.6-0.8), excellent (>0.8) and prefect agreement (+1) [60,61].

Prediction accuracy of models and variable importance
The score of individual models based on average AUC ranged from 0.86 (MAXENT.Phillips) to 0.97 (GBM/RF); while the average TSS scores ranged from 0.73 (MAXENT. Phillips) to 0.92 (GBM) ( Table 1). Based on the AUC score, the best performing model was the GBM and the RF while based on the TSS score, it was GBM. Nonetheless, all the single models performed better than random models but not better than an ensemble model ( Figure 2). The median AUC of the ensemble model was 0.99 while TSS was 0.97. As predicted by the model, compared to the current period, the suitable habitat area of parthenium weed will increase by the years 2050 and 2070. Amongst the nine predictor variables used for modelling, mean temperature of the wettest quarter (Bhutan_Bio8) was the most important variable influencing parthenium weed distribution in Bhutan. Precipitation seasonality (Bhu-tan_Bio15), precipitation of the warmest quarter (Bhu-tan_Bio18) and the min temperature of the coldest month (Bhutan_Bio6) influenced this to some extent ( Table 2). Parthenium weed suitability increased with increasing mean temperature of wettest quarter (Bhutan_Bio8) and precipitation seasonality (Bhutan_Bio15) while it decreased with increasing precipitation of warmest quarter (Bhutan_Bio18) and the min temperature of the coldest month (Bhutan_Bio6) (Figure 3).

Current and future distribution of parthenium weed in Bhutan
The potential distribution of parthenium weed in Bhutan based on the current climatic conditions and occurrences points showed that about 1,099.1 km 2 (2.83%) of the total area of Bhutan is suitable for parthenium weed invasion ( Figure 4 and  (Figure 4 and Table 4). In terms of loss, the highest loss in suitability was determined under RCP4.5 2070 with about 178.35 km 2 while the least was determined under RCP4.5 2050 ( Figure 4 and Table 4). Similarly, the highest stability was predicted under RCP4.5 2050 with about 1,130.31 km 2 while the lowest was predicted under RCP4.5 2070 with about 957.53 km 2 ( Figure 4 and Table 4). The difference in suitability in terms of elevation range at different climate periods was also determined. There was a stark difference in the elevation range of parthenium weed distribution between current and future climates. The occurrence points were collected within the elevation range between 210-1,423 m asl but ensemble model predicted the elevation range for current distribution to be between 109-2,178 m asl. The upper elevation range increased under both the future climates and scenarios and ranged between 2,513-2,931 m asl. The highest elevation of about 2,931 m asl was predicted under RCP8.5 2070; an upward migration by about 753 m asl as compared to the current distribution.

Discussion
While species distribution modelling is a popular choice to assess potential risk of species invasion [62], there are various factors that influence the performance and thereby the outcome of the models. These include the extent and resolution of the study area and threshold used [63,64], number of occurrences used for model training [65], species type -whether common or rare [66], variable selection technique, environmental predictors and species distributional data [67], choice of models [68], modelling approaches [69], GCMs used [70], model complexity [71], collinearity [72], presence-only or presence/ absence data used [62,73], method and number of back-ground points used [55,62], method of evaluation [74] and binary maps' thresholds used [75]. Several algorithms are recommended for modelling as no single model performs best under all situations [76]. Also, ensemble models account for uncertainties in predictions of different algorithms and performs better than a single modelling algorithm [52,36]. Therefore, we decided to use ensemble modelling approach in our study. We used both the AUC and TSS to evaluate models as the use of single AUC to pick a model is not recommended [60]. The performance of individual model ranged from good to excellent based on the average AUC, which ranged from 0.86 to 0.97. True skills statistics ranged from 0.73 to 0.92 indicating useful to excellent performance. The prediction accuracy was significantly improved by combining models and building an ensemble in which the AUC was 0.99 (excellent) and TSS was 0.97 (excellent). Given that Bhutan values its biodiversity conservation and environmental protection, ensuring both are mainstreamed into the country's development programs and policies, the use of species distribution modelling seems inevitable. Yet, species distribution modelling has never been used, especially in spatial planning and deci-sion-making in the context of plant invasion [12] tried to model six invasive alien species, including parthenium weed in Bhutan using the most commonly implemented MaxEnt tool [34]. However, the study did not provide appropriate details in terms of districts suitability and elevational changes other than range expansion, contraction, and directional change, which have limited use when planning for management programs. Similarly, using the MaxEnt tool, [33] studied parthenium weed distribution and its impact in the West-Central region but did not provide much detail. Also, they considered only one future period (2070). We have modelled for the first time the distribution of parthenium weed in Bhutan based on ensemble modelling approach using the BIOMOD2 package developed by [35] and [36] in R v 4.0.3 [37] for the current and two future climate periods (2050 and 2070) ( Figure 5). In addition to providing the details of invasion at the district level, we have also analyzed the elevation migration of parthenium weed to better determine the impact of climate change in this mountainous country. Our study will help direct future management plans and inform policy makers on the need for early detection and management in newly invaded areas. Besides climatic variables, there are other factors such as dispersal, demography and species interactions that determine species distribution [77]. However, information on such factors is not easily available and including them increases the complexity of the modelling process [78]. Additionally, in mountainous country like Bhutan where elevation ranges from 150 m asl to more than 7,000 m asl with steep slopes, topography can play a very critical role in creating barriers to the dispersal of species. Therefore, we included slope as one of the predictor variables in the modelling process. However, we did not find any influence of slope in parthenium weed distribution in Bhutan, which did not agree with the findings of [23] in Nepal, where slope was found to have contributed the most. This was attributed to the choice of model algorithm where we have used an ensemble modelling approach instead of the single model modelling approach. If individual model was considered, MaxEnt did give slope as the most influencing variable ( Table 2). Climate change has expanded the ranges of species [31,79]. Like other studies in the region [31,44,80,81], we have also demonstrated the impact of climate change and showed that the land suitability for parthenium weed invasion significantly increased under future climatic conditions. The percentage increase ranged from 7.39 to 13.95% of Bhutan's total area, or a 161 to 393% increase of the current predicted potential suitability. These findings are in sharp contrast to those of [12] where they have reported suitable land for parthenium weed invasion to contract in 2070 by about 4.15% as compared to the present 8.36%. Again, this may be attributed to the modelling method and choice of algorithm used, as they only used the MaxEnt program. Under current climatic conditions, out of 20 districts, 17 districts showed suitability to parthenium weed invasion. In the future, because of climate change, two additional districts (Thimphu and Paro) become suitable. Bumthang remained free from invasion under all the climatic conditions tested. These findings were also validated with the present day occurrence points and with the observations of agriculture extension officials working in those districts. Occurrence points were collected only from the western and southern districts of Bhutan but the prediction showed invasion suitability in eastern and central districts where occurrence data were not collected such as Lhuentse, Mongar, Trashigang, Trongsa and Zhemgang. Additionally, [28]  Among the nine predictor variables used in the modelling, the most important ones that explained the species' environmental requirements best were the mean tempera-ture of the wettest quarter (Bhutan_Bio8), precipitation seasonality (Bhutan_Bio15), precipitation within the warmest quarter (Bhutan_Bio18) and the min temperature of the coldest month (Bhutan_Bio6). The suitability of parthenium weed increased with increasing mean temperature of the wettest quarter (Bhutan_Bio8) and with precipitation seasonality (Bhutan_Bio15) while it decreased with increasing precipitation of the warmest quarter (Bhutan_Bio18) and the min temperature of the coldest month (Bhutan_Bio6). The low temperature of the mountains is the main limiting factor preventing species dispersal from lowlands into the mountains [82,83]. According to the observations of [84] and [85] in Australia, the distribution of parthenium weed was limited to the temperatures between 5°C and 40°C. One of the major determinants of the parthenium weed distribution is moisture availability [86,87]. However, climate change will create climatically suitable niches in higher elevations through increased temperature and parthenium weed will thrive as a result, and as predicted by our model. Climate projections by the National Center for Hydrology and Meteorology (NCHM) in Bhutan have predicted an increased surface temperature of between 0.8 to 1.6°C under RCP4.5 and between 0.8 to >3.2°C under RCP 8.5 during the period 2021 to 2050 [42]. We have predicted a maximum upward elevation shift to 2,931 m asl; about 753 m asl increase under RCP8.5 and by 2070. Presence of parthenium weed at 2,600 m asl was reported in Bhutan [13], while in Ethiopia, the highest altitude record of infestation was at 2,627 m [88]. The upslope movement of parthenium weed would mean highland pastures coming under threat from invasion and making such pastures unavailable for livestock production, especially for high altitude dwelling animals like yaks. Besides climate change, anthropogenic disturbance also drive the plant invasion in mountain ecosystems [89,83]. Bhutan has undergone unprecedented socio-economic development since its transition to democracy in 2008. Road network expansion is one of the main pledges of the political parties. Even those areas where construction was once restricted (e.g., Sakteng Wildlife Sanctuary) is now connected by roads. As roads are considered to be the major conduits through which parthenium weed is dispersed [90,88], roads will invariably allow parthenium weed to invade not only in lowland regions but also at higher elevations in Bhutan. Anthropogenic disturbance creating suitable habitat for parthenium weed may already be occurring as seen from the location of occurrence points collected. The highest infestation and occurrences were collected along the roads where one of the biggest hydropower projects in the country, Punatshangchu Hydro-Power Project II is currently being built. It may be conjectured here that excavations and extensive movement of vehicles carrying hydro-electricity equipment might have facilitated the spread of parthenium weed. It is more challenging to control invasive weeds in the mountainous areas than in adjoining lowland areas [91]. Therefore, early detection and the prevention of parthenium weed population establishment will be very critical. Bhutan is renowned for its very strong environmental protection policy, high density of plants and animals and the constitution necessitates that the country maintain 60% of the land under forest cover in perpetuity [92]. Given the unfettered support for environmental protection, the country's mountainous terrain with steep slopes, high elevations, and dense forest may all contribute to preventing encroachment of invasive species. However, as stated, anthropogenic disturbances are increasing as political parties try to fulfill their pledges, and as the country strives to alleviate poverty and secure food for all. Further, Bhutan's dependence on agriculture for economy and livelihood while only about 7% of the country's area is arable, biological invasion may have a significant impact on food security and food production. Farmers have already reported impact of biological invasions in Bhutan [20]. Thirty horses have succumbed to Crofton weed consumption [19]. In Royal Manas National Park, 61 ha of forest that was cleared to convert to grassland for wildlife has been invaded by Crofton weed [14]. Furthermore, Bhutan is highly vulnerable to climate change [42]. While in-depth studies on the impacts of IAPS in Bhutan has not been undertaken, inventories are being initiated [10,93,16]. In this present study we have modelled the distribution of one of the most invasive weeds, parthenium weed in Bhutan under the current and future climates and have shown that it will become highly invasive as a result of climate change and with the species even moving into higher elevation zones. This study should encourage land managers into undertaking early detection, monitoring and developing management plans for parthenium weed in Bhutan. One of the management options that could be explored is the use of biological agents, an approach that is in line with Bhutan's policy of restricted chemical use for environmental protection. The presences of two biological agents, the Mexican beetle (Zygogramma bicolorata Pallister, 1953) and the winter rust (Puccinia abrupta var. partheniicola (H.S.Jacks.) Parmelee) have already arrived in Bhutan [94]. These agents have been introduced and have successfully established in Australia [95] and India [96] and are an important component of those countries management plans of parthenium weed.

Conclusion
We have modelled the impact of climate change on the distribution of one of the most invasive weeds, parthenium weed in Bhutan. Through our study, we have shown that under the future climate period (2050 and 2070) the suitable habitat areas for parthenium weed will expand by as much as five times the current distribution, and come to infest all districts of Bhutan, except for the high mountainous district of Bumthang. Under the future climate, the highest land suitability for parthenium weed is predicted under RCP4.5 2050 period with about 5,419.69 km 2 . Generally, the districts in west and south showed more suitability than the districts in the east and central. The highest elevation of suitability was predicted at 2,931 m asl; an upward shift of about 753 m by 2070. These observations suggest that many of the country's native species will be at risk from displacement by parthenium weed and might even be pushed to extinction. Food security and food production in Bhutan will also be constrained as parthenium weed commonly invades agricultural areas. Therefore, early detection and eradication plan of new parthenium weed populations will be very important as currently, only about 2.83% of the country's area is predicted to have suitability to parthenium weed invasion, which appears to be quite manageable. Early detection and rapid response should remain our priority since delay in control programs will be very costly with eradication impossible in future.