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

Sangay Dorji, Department of Agriculture and Food Sciences, The University of Queensland, Brisbane, Australia, Email: s.dorji@uqconnect.edu.au

Received: 15-Mar-2022, Manuscript No. JENVOH-22-57196; Editor assigned: 17-Mar-2022, Pre QC No. JENVOH-22-57196 (PQ); Reviewed: 01-Apr-2022, QC No. JENVOH-22-57196; Revised: 06-Apr-2022, Manuscript No. JENVOH-22-57196 (R); Published: 13-May-2022

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², 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-16], including Crofton weed (Ageretina adenophora (Spreng.) King & H.Rob.), Siam weeds (Chromolaena odorata L.) R.M.King & H.Rob.), lantana (Lantana camara L.), mile-a-minute (Mikania micrantha Kunth) and parthenium weed (Parthenium hysterophorus L.). 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 BIOMOD2 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.

Materials and Methods

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.

Figure 1. Occurrence points (red dots) along the surveyed roads (green). Note: ( ) Occurrence, ( ) surveyed roads, ( ) Highway ( ) Districts

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²) 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-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: (i) Mean Diurnal Range (Bhutan_Bio2), (ii) Isothermality (Bhutan_Bio3), (iii) Min Temperature of Coldest Month (Bhutan_Bio6), (iv) Mean Temperature of Wettest Quarter (Bhutan_Bio8), (v) Precipitation Seasonality (Bhutan_ Bio15), (vi) Precipitation of Warmest Quarter (Bhutan_ Bio18), (vii) Precipitation of the Coldest Quarter (Bhutan_Bio19), (viii) Aspect, and (ix) Slope. 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 (MAXENT. 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 StaStatistics (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].

Results

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 (Bhutan_ Bio15), precipitation of the warmest quarter (Bhutan_ 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).

Figure 2. Performance of individual model and ensemble model based on the Receiver Operating Characteristics (ROC) and the Area Under Curve (AUC) and True Skills Statistics (TSS) in predicting the distribution of P. hysterophorus in Bhutan. Models used: Generalized Additive Model (GAM), Generalized Linear Model (GLM), Generalized Boosted Model (GBM), Random Forest (RF), and Maximum Entropy (MAXENT.Phillips). The bottom of the whisker of the box-plot represents minimum, bottom edge represents lower quartile, bold line inside represents median, upper edge represents upper quartile, top of the whisker represents maximum and filled dots represents outliers. The predictive performance of the ensemble model (AUC = 0.99, TSS = 0.97) is represented by the coloured dashed line.

Figure 3. Response curve of probability of occurrence obtained for top four contributing predictor variables of ensemble model by weighted-mean approach based on Area Under Curve (AUC). Precipitation seasonality (Bhutan_Bio15), precipitation of warmest quarter (Bhutan_Bio18), min temperature of the coldest month (Bhutan_Bio6), and mean temperature of wettest quarter (Bhutan_Bio8).

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.83%) of the total area of Bhutan is suitable for parthenium weed invasion (Figure 4 and Table 3). This includes regions within 17 out of the 20 districts of Bhutan. Amongst districts, Wangduephodrang had the highest suitability of about 302.09 km2 (27.49%) followed by Punakha, Dagana and Samtse with 142.68 km² (12.98%), 129.31 km² (11.76%) and 107.01 km² (9.74%), respectively. Gasa and Haa had the lowest suitability.

Figure 4. Ensemble maps showing current and future distribution of P. hysterophorus in Bhutan. Leaf green areas indicate suitable areas while white areas indicate not suitable areas. The maps were prepared in ArcGIS 10.4. Note: ( ) Suitable, ( ) Not suitable, ( ) Districts

Under future climatic projections, the land suitable for parthenium weed invasion increased in both the climate periods 2050 and 2070 (Figure 4 and Table 3). The highest suitability was predicted under RCP4.5 2050 with about 5,419.69 km² (13.95% of Bhutan’s area), followed by RCP8.5 2070 with about 3,954.97 km² (13.01% of Bhutan’s area). The lowest is predicted under RCP4.5 2070 with about 1,771.27 km² (Table 3). There was also a change in the number of districts suitable to parthenium weed invasion. The 17 districts that were suitable under current remained stable under RCP4.5 2050 scenario while additional two districts (Paro and Thimphu) became suitable under RCP8.5 2050. However, under RCP4.5 2070 scenario, only one additional district (Paro) became suitable compared to current prediction while under RCP8.5 2070, districts’ suitability remained same as under RCP8.5 2050 scenario. Bumthang remained free from parthenium weed invasion under any period.

District-wise, under RCP4.5 2050 scenario, Sarpang had the highest suitability with about 706.72 km² followed by Dagana with about 692.23 km2. Under RCP8.5 2050 scenario, Dagana had the highest suitability with about 462.60 km² followed by Wangduephodrang with about 429 km². Under RCP4.5 2070 scenario, Dagana had the highest suitability with about 441.42 km² followed by Wangduephodrang with about 434.74 km². Under RCP8.5 2070 scenario, the highest suitability was predicted in Dagana with about 681.08 km² followed by Wangduephodrang with about 560.69 km².

In terms of change in suitability between current and future climate periods, the highest gain in suitability was determined under RCP4.5 2050 scenario with about 4,391.94 km² while the lowest was determined under RCP 4.5 2070 scenario with about 1,959.65 km² (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² 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² while the lowest was predicted under RCP4.5 2070 with about 957.53 km² (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 decision-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.

Figure 5. Change in the distribution of P. hysterophorus between current and future climate scenarios. The maps were prepared in ArcGIS 10.4. Note: ( ) Gain, ( ) loss, ( ) stable, ( ) Not suitable, ( ) Districts.

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] reported parthenium weed in Mongar, Trashigang and Trongsa while the first author detected five plants of parthenium weed in Thimphu around his residence area (27.51396°N, 89.6404°E, 2, 400 m asl) (model predicted suitability only in RCP8.5 2050 and 2070).

Among the nine predictor variables used in the modelling, the most important ones that explained the species’ environmental requirements best were the mean temperature 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 km2. 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.

Acknowledgements

We would like to thank Mr. Sangay Wangdi, Department of Agriculture, Ministry of Agriculture and Forests, Thimphu, Bhutan for sharing local names and his experience with the weed.

Authors’ Contributions

Sangay Dorji: Conceptualization, Survey Design, Field Work, Data Collection, Analysis, Writing–original draft. Lakey Lakey: Writing–original draft, including review and editing. Tshering Wangchen: Writing–original draft, including review and editing. Steve Adkins: Supervision, Manuscript Review and Editing.

Competing Interests

The authors declare that they have no competing interests.

References

1. Hooper DU, Chapin Iii FS, Ewel JJ, Hector A, Inchausti P, Lavorel S, et al. Effects of biodiversity on ecosystem functioning: a consensus of current knowledge. Ecol Monogr 2005; 75:3-35.

   [Crossref] [Google scholar]

2. Blackburn TM, Bellard C, Ricciardi A. Alien versus native species as drivers of recent extinctions. Front Ecol Environ 2019; 17:203-207.

   [Crossref] [Google scholar]

3. Downey PO, Richardson DM. Alien plant invasions and native plant extinctions: a six-threshold framework. AoB Plants 2016; 8:plw047.

   [Crossref] [Google scholar] [Pubmed]

4. Bongaarts J. IPBES, 2019. Summary for policymakers of the global assessment report on biodiversity and ecosystem services of the Intergovernmental Science‐Policy Platform on Biodiversity and Ecosystem Services. Popul Dev Rev 2019; 45:680-681.

   [Crossref] [Google scholar]

5. Finch DM, Butler JL, Runyon JB, Fettig CJ, Kilkenny FF, Jose S, et al. Effects of Climate Change on Invasive Species. In:  Poland TM, Patel-Weynand T, Finch DM, Miniat CF, Hayes DC and Lopez VM  (ed) Invasive Species in Forests and Rangelands of the United States: A Comprehensive Science Synthesis for the United States Forest Sector. Springer International Publishing, Cham. 2021; 57-83.

   [Crossref] [Google scholar]

6. Hellmann JJ, Byers JE, Bierwagen BG, Dukes JS. Five Potential Consequences of Climate Change for Invasive Species. Conserv Biol 2008; 22:534-543.

   [Crossref] [Google scholar] [Pubmed]

7. Effects of biodiversity on ecosystem functioning: a consensus of current knowledge
8. Rai K, Singh JS. Invasive alien plant species: Their impact on environment, ecosystem services and human health. Ecol Indic 2020; 111:106020.

   [Crossref] [Google scholar] [Pubmed]

9. Crossref
10. Dorjee, Johnson SB, Buckmaster AJ, Downey PO. Weeds in the land of Gross National Happiness: Knowing what to manage by creating a baseline alien plant inventory for Bhutan. Biol Invasions 2020; 22:2899-2914.

    [Crossref] [Google scholar]

11. Mahat K, Mitchell A, Zangpo T. An updated global COI barcode reference data set for Fall Armyworm (Spodoptera frugiperda) and first record of this species in Bhutan. J Asia-Pac Entomol 2021; 24:105-109.

    [Crossref] [Google scholar]

12. Thinley U, Banterng P, Gonkhamdee S, Katawatin R. Distributions of Alien Invasive Weeds under Climate Change Scenarios in Mountainous Bhutan. Agronomy 2019; 9:442.

    [Crossref] [Google scholar]

13. Tshering C, Adkins S. Parthenium weeds in Bhutan. In: International Parthenium News, 2012; 5:9-10.
14. Wangdi T. Grassland in Manas park under weed attack. Kuensel 2015.
15. Weiss J, Thinley T, Nidup K, Ghimiray M, Wandi S, Dochen T. Working with weeds in the Land of the Thunderdragon: an opportunity to prevent weed invasions. In: Proceedings of the 14th Australian Weeds Conference. 2004.

    [Google scholar]

16. Yangzom R, Dorji K, Dorji T, Dorji R, Wangmo C, Gyeltshen C. A Pictorial Guide to Major Invasive Plant Species of Bhutan. National Biodiversity Center 2018.
17. NBC (2014) National Biodiversity Strategies and Action Plan of Bhutan. National Biodiversity Center. 2014.

    [Google scholar]

18. Kaur M, Aggarwal NK, Kumar V, Dhiman R. Effects and Management of Parthenium hysterophorus: A Weed of Global Significance. Int Sch Res Notices 2014:368647.

    [Crossref] [Google scholar] [Pubmed]

19. Dema T. Equine fatalities from epotorium plant. In: Kuenselonline 2009.

    [Google scholar]

20. Suberi B, Tiwari K, Gurung D, Bajracharya R, Sitaula B. People’s perception of climate change impacts and their adaptation practices in Khotokha valley, Wangdue, Bhutan. Indian J Tradit Know 2018; 17:97-105.

    [Google scholar]

21. Pallewatta N, Reaser J, Gutiérrez A. Invasive alien species in South-Southeast Asia: national reports and directory of resources. Global Invasive Species Programme 2003.

    [Google scholar]

22. Hoffmann BD, Broadhurst LM. The economic cost of managing invasive species in Australia. NeoBiota 2016.  

    [Crossref] [Google scholar]

23. Maharjan S, Shrestha BB, Joshi MD, Devkota A, Muniappan R, Adiga A et al. Predicting suitable habitat of an invasive weed Parthenium hysterophorus under future climate scenarios in Chitwan Annapurna Landscape, Nepal. J Mt Sci 2019; 16:2243-2256.

    [Crossref] [Google scholar]

24. Hiremath AJ, Krishnan S. India Knows Its Invasive Species Problem But This Is Why Nobody Can Deal With it Properly. The Wire 2016.

    [Google scholar]

25. Adkins S, McClay A, Bajwa AA. Biology and Ecology. CAB International 2019; 7-39.

    [Crossref] [Google scholar]

26. Bajwa AA, Chauhan BS, Farooq M, Shabbir A, Adkins SW. What do we really know about alien plant invasion? A review of the invasion mechanism of one of the world’s worst weeds. Planta 2016; 244:39-57.

    [Crossref] [Google scholar] [Pubmed]

27. Shabbir A, McConnachie A, Adkins SW. Parthenium weed : biology, ecology and management. AGRIS 2019.

    [Google scholar]

28. Parker C. Weeds of Bhutan. National Plant Protection Center 1992.

    [Google scholar]

29. Shrestha BB, Shabbir A, Adkins SW. Parthenium hysterophorus in Nepal: a review of its weed status and possibilities for management. Weed Res 2015; 55:132-144.

    [Crossref] [Google scholar]

30. Nguyen T, Bajwa AA, Navie S, O’Donnell C, Adkins S. Parthenium weed (Parthenium hysterophorus L.) and climate change: the effect of CO2 concentration, temperature, and water deficit on growth and reproduction of two biotypes. Environ Sci and Pollut R 2017; 24:10727-10739.

    [Crossref] [Google scholar] [Pubmed]

31. Bellard C, Thuiller W, Leroy B, Genovesi P, Bakkenes M, Courchamp F. Will climate change promote future invasions? Glob Change Biol 2013; 19:3740-3748.

    [Crossref] [Google scholar] [Pubmed]

32. Thinley U, Banterng P, Gonkhamdee S, Katawatin R. Distributions of Alien Invasive Weeds under Climate Change Scenarios in Mountainous Bhutan. Agronomy 2019; 9:442.

    [Crossref] [Google scholar]

33. Chhogyel N, Kumar L, Bajgai Y. Invasion status and impacts of parthenium weed (Parthenium hysterophorus) in West-Central region of Bhutan. Biol Invasions 2021.

    [Crossref] [Google scholar]

34. Santini L, Benítez-López A, Maiorano L, Čengić M, Huijbregts MAJ. Assessing the reliability of species distribution projections in climate change research. Divers Distrib 2021; 27:1035-1050.

    [Crossref] [Google scholar]

35. Thuiller W, Georges D, Engler R, Breiner F, Georges M, Thuiller C. Package ‘biomod2’. 2016.

    [Google scholar]

36. Thuiller W, Lafourcade B, Engler R, Araújo MB. BIOMOD–a platform for ensemble forecasting of species distributions. Ecography 2009; 32:369-373.

    [Crossref] [Google scholar]

37. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing. 2020.

    [Google scholar]

38. Brown JL, Bennett JR, French CM. SDMtoolbox 2.0: the next generation Python-based GIS toolkit for landscape genetic, biogeographic and species distribution model analyses. PeerJ 2017; 5:e4095.

    [Crossref] [Google scholar] [Pubmed]

39. Fick SE, Hijmans RJ. WorldClim 2: new 1‐km spatial resolution climate surfaces for global land areas. Int J Climatol 2017; 37:4302-4315.

    [Crossref] [Google scholar]

40. Karger DN, Conrad O, Böhner J, Kawohl T, Kreft H, Soria-Auza RW, et al. Climatologies at high resolution for the earth's land surface areas. Sci Data 2017; 4:170122.

    [Crossref] [Google scholar] [Pubmed]

41. Google scholar
42. Dorji T, Tamang TB. Report on the analysis of historical climate and climate projection for Bhutan. National Center for Hydrology and Meteorology, Royal Government of Bhutan, PO Box: 2017, Thimphu, Bhutan. NCHM 2019; 1-52.

    [Google scholar]

43. Aryal A, Shrestha UB, Ji W, Ale SB, Shrestha S, Ingty T, et al. Predicting the distributions of predator (snow leopard) and prey (blue sheep) under climate change in the Himalaya. Ecol Evol 2016; 6:4065-4075.

    [Crossref] [Google scholar] [Pubmed]

44. Lamsal P, Kumar L, Aryal A, Atreya K. Invasive alien plant species dynamics in the Himalayan region under climate change. Ambio 2018; 47:697-710.

    [Crossref] [Google scholar] [Pubmed]

45. Mishra V, Kumar D, Ganguly AR, Sanjay J, Mujumdar M, Krishnan R, et al. Reliability of regional and global climate models to simulate precipitation extremes over India. JGR Atmospheres 2014; 119:9301-9323.

    [Crossref] [Google scholar]

46. Sharmila S, Joseph S, Sahai AK, Abhilash S, Chattopadhyay R. Future projection of Indian summer monsoon variability under climate change scenario: An assessment from CMIP5 climate models. Glob Planet Change 2015; 124:62-78.

    [Crossref] [Google scholar]

47. Su J, Aryal A, Nan Z, Ji W. Climate Change-Induced Range Expansion of a Subterranean Rodent: Implications for Rangeland Management in Qinghai-Tibetan Plateau. PLoS ONE 2015; 10:e0138969.

    [Crossref] [Google scholar] [Pubmed]

48. Merow C, Smith MJ, Silander Jr JA. A practical guide to MaxEnt for modeling species' distributions: what it does, and why inputs and settings matter. Ecography 2013; 36:1058-1069.

    [Crossref] [Google scholar]

49. Naimi B, Hamm NAS, Groen TA, Skidmore AK, Toxopeus AG. Where is positional uncertainty a problem for species distribution modelling? Ecography 2014; 37:191-203.

    [Crossref] [Google scholar]

50. Menard SW. Applied logistic regression analysis. Sage Publications London. 2001.

    [Google scholar]

51. Elith J, Leathwick JR. Species Distribution Models: Ecological Explanation and Prediction Across Space and Time. Annu Rev Ecol Evol Syst 2009; 40:677-679.

    [Crossref] [Google scholar]

52. Aguirre-Gutiérrez J, van Treuren R, Hoekstra R, van Hintum TJL. Crop wild relatives range shifts and conservation in Europe under climate change. Divers Distrib 2017; 23:739-750.

    [Crossref] [Google scholar]

53. Becker EA, Carretta JV, Forney KA, Barlow J, Brodie S, Hoopes R, et al. Performance evaluation of cetacean species distribution models developed using generalized additive models and boosted regression trees. Ecol Evol 2020; 10:5759-5784.

    [Crossref] [Google scholar] [Pubmed]

54. Elith J, Graham CH, P. Anderson R, Dudík M, Ferrier S, Guisan A, et al. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 2006;29:129–151.

    [Crossref] [Google scholar]

55. Barbet‐Massin M, Jiguet F, Albert CH, Thuiller W. Selecting pseudo‐absences for species distribution models: how, where and how many? Methods Ecol Evol 2012; 3:327-338.

    [Crossref] [Google scholar]

56. Araújo MB, Pearson RG, Thuiller W, Erhard M. Validation of species-climate impact models under climate change. Glob Change Biol 2005;11:1504-1513.

    [Crossref] [Google scholar]

57. Fielding AH, Bell JF. A review of methods for the assessment of prediction errors in conservation presence/absence models. Envir Conserv 1997; 24:38-49.

    [Crossref] [Google scholar]

58. Peterson A, Soberón J, Pearson R, Anderson R, Martínez-Meyer E, Nakamura M, et al. Ecological Niches and Geographic Distributions. Princeton University Press 2011.

    [Crossref] [Google scholar]

59. Swets JA. Measuring the accuracy of diagnostic systems. Science 1988; 240:1285-1293.

    [Crossref] [Google scholar] [Pubmed]

60. Allouche O, Tsoar A, Kadmon R. Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). J Appl Ecol 2006;43:1223-1232.

    [Crossref] [Google scholar]

61. Beaumont LJ, Graham E, Duursma DE, Wilson PD, Cabrelli A, Baumgartner JB. Which species distribution models are more (or less) likely to project broad-scale, climate-induced shifts in species ranges? Ecol Modell 2016; 342:135-146.

    [Crossref] [Google scholar]

62. Elith J. Predicting Distributions of Invasive Species. In: Robinson AP, Walshe T, Burgman MA and Nunn M (ed) Invasive Species Risk Assessment and Management. Cambridge University Press, London. 2017; 93-129.

    [Google scholar]

63. Acevedo P, Jiménez-Valverde A, Lobo JM, Real R. Delimiting the geographical background in species distribution modelling. J Biogeogr 2012; 39:1383-1390.

    [Crossref] [Google scholar]

64. Vale CG, Tarroso P, Brito JC. Predicting species distribution at range margins: testing the effects of study area extent, resolution and threshold selection in the Sahara–Sahel transition zone. Divers Distrib  2014; 20:20-33.

    [Crossref] [Google scholar]

65. Guisan A, Graham CH, Elith J, Huettmann F. Sensitivity of predictive species distribution models to change in grain size. Divers Distrib 2007; 13:332-340.

    [Crossref] [Google scholar]

66. Norberg A, Abrego N, Blanchet FG, Adler FR, Anderson BJ, Anttila J, et al. A comprehensive evaluation of predictive performance of 33 species distribution models at species and community levels. Ecol Monogr 2019; 89:e01370.

    [Crossref] [Google scholar]

67. Stankowski PA, Parker WH. Species distribution modelling: Does one size fit all? A phytogeographic analysis of Salix in Ontario. Ecol Modell 2010; 221:1655-1664.

    [Crossref] [Google scholar]

68. Benito BM, Cayuela L, Albuquerque FS, O'Hara RB. The impact of modelling choices in the predictive performance of richness maps derived from species‐distribution models: guidelines to build better diversity models. Methods Ecol Evol 2013; 4:327-335.

    [Crossref] [Google scholar]

69. Kaky E, Nolan V, Alatawi A, Gilbert F. A comparison between Ensemble and MaxEnt species distribution modelling approaches for conservation: A case study with Egyptian medicinal plants. Ecol Inform 2020; 60:101150.

    [Crossref] [Google scholar]

70. Steen V, Sofaer H, Skagen SK, Ray AJ, Noon BR. Projecting species’ vulnerability to climate change: Which uncertainty sources matter most and extrapolate best? Ecol Evol 2017;7:8841-8851.

    [Crossref] [Google scholar] [Pubmed]

71. Brun P, Thuiller W, Chauvier Y, Pellissier L, Wüest RO, Wang Z, et al. Model complexity affects species distribution projections under climate change. J Biogeogr 2020; 47:130-142.

    [Crossref] [Google scholar]

72. Dormann CF. Effects of incorporating spatial autocorrelation into the analysis of species distribution data. Glob Ecol Biogeogr 2007; 16:129-138.

    [Crossref] [Google scholar]

73. Glor RE, Warren D. Testing ecological explanations for biogeographic boundaries. Evolution 2011; 65:673-683.

    [Crossref] [Google scholar] [Pubmed]

74. Lobo JM, Jiménez-Valverde A, Real R. AUC: a misleading measure of the performance of predictive distribution models. Glob Ecol Biogeogr 2008; 17:145-151

    [Crossref] [Google scholar]

75. Liu C, White M, Newell G. Selecting thresholds for the prediction of species occurrence with presence-only data. J Biogeogr 2013; 40:778-789.

    [Crossref] [Google scholar]

76. Qiao H, Soberón J, Peterson AT, Kriticos D. No silver bullets in correlative ecological niche modelling: insights from testing among many potential algorithms for niche estimation. Methods Ecol Evol 2015; 6:1126-1136. [Crossref]

    [Google scholar]

77. Zurell D. Integrating demography, dispersal and interspecific interactions into bird distribution models. J Avian Biol 2017; 48:1505-1516.

    [Crossref] [Google scholar]

78. Evans TG, Diamond SE, Kelly MW. Mechanistic species distribution modelling as a link between physiology and conservation. Conserv Physiol 2015; 3:cov056-cov056.

    [Crossref] [Google scholar] [Pubmed]

79. Taylor S, Kumar L, Reid N, Kriticos DJ. Climate Change and the Potential Distribution of an Invasive Shrub, Lantana camara L. PLoS ONE 2012; 7:e35565.

    [Crossref] [Google scholar] [Pubmed]

80. Mainali KP, Warren DL, Dhileepan K, McConnachie A, Strathie L, Hassan G, et al. Projecting future expansion of invasive species: comparing and improving methodologies for species distribution modeling. Glob Change Biol 2015; 21:4464-4480.

    [Crossref] [Google scholar] [Pubmed]

81. Thapa S, Chitale V, Rijal SJ, Bisht N, Shrestha BB. Understanding the dynamics in distribution of invasive alien plant species under predicted climate change in Western Himalaya. PLoS ONE 2018; 13: e0195752.

    [Crossref] [Google scholar] [Pubmed]

82. Alexander JM, Kueffer C, Daehler CC, Edwards PJ, Pauchard A, Seipel T. Assembly of nonnative floras along elevational gradients explained by directional ecological filtering. Proc Natl Acad Sci USA 2011; 108:656.

    [Crossref] [Google scholar] [Pubmed]

83. Petitpierre B, McDougall K, Seipel T, Broennimann O, Guisan A, Kueffer C. Will climate change increase the risk of plant invasions into mountains? Ecol Appl 2016; 26:530-544.

    [Crossref] [Google scholar] [Pubmed]

84. Dale IJ. Parthenium weed in the Americas: a report on the ecology of Parthenium hysterophorus in South, Central and North America. Australian Weeds 1981; 1:8-14.

    [Google scholar]

85. Doley D. Parthenium weed (Parthenium hysterophorus L.): Gas exchange characteristics as a basis for prediction of its geographical distribution. Aust J Agric Res 1977; 28:449-460.

    [Crossref] [Google scholar]

86. Ayele S. The Impact of Parthenium (Parthenium hysterophorus L.) on the Range Ecosystem Dynamics of the Jijiga Rangeland, Ethiopia. 2007.

    [Google scholar]

87. Goodall J, Braack M, de Klerk J, Keen C. Study on the early effects of several weed-control methods on Parthenium hysterophorus L. Afr J Range Forage Sci 2010;27: 95-99.

    [Crossref] [Google scholar]

88. McConnachie AJ, Strathie LW, Mersie W, Gebrehiwot L, Zewdie K, Abdurehim A, et al. Current and potential geographical distribution of the invasive plant Parthenium hysterophorus (Asteraceae) in eastern and southern Africa. Weed Res 2011; 51:71-84.

    [Crossref] [Google scholar]

89. Davis MA, Grime JP, Thompson K. Fluctuating Resources in Plant Communities: A General Theory of Invasibility. J Ecol 2000; 88:528-534.

    [Crossref] [Google scholar]

90. Blackmore P, Johnson S. Continuing successful eradication of parthenium weed (Parthenium hysterophorus) from New South Wales, Australia. In: Seventeenth Australasian Weeds Conference, Christchurch, New Zealand. NSW 2010; 2:382-385.

    [Google scholar]

91. Pauchard A, Milbau A, Albihn A, Alexander J, Burgess T, Daehler C, et al. Non-native and native organisms moving into high elevation and high latitude ecosystems in an era of climate change: new challenges for ecology and conservation. Biol Invasions 2016; 18:345-353.

    [Crossref] [Google scholar]

92. RGoB. The Constitution of the Kingdom of Bhutan. 2008.

    [Google scholar]

93. Gyeltshen C, Dema S, Northey D, Prasad K. Biodiversity Statistics of Bhutan 2017 - A Preliminary Baseline. National Biodiversity Center, Thimphu 2019.

    [Google scholar]

94. Dorji S, Adkins S. First record of the beetle Zygogramma bicolorata and the rust fungus Puccinia abrupta var. partheniicola on Parthenium hysterophorus in Bhutan. International Parthenium News. 2020; 15:1-17.
95. Dhileepan K. Biological Control of Parthenium (Parthenium hysterophorus) in Australian Rangeland Translates to Improved Grass Production. Weed Sci 2007; 55:497-501.

    [Crossref] [Google scholar]

96. Kumar S. Biological control of Parthenium in India: Status and prospects. Indian J Weed Sci (Online) 2009; 41:1-18.

    [Crossref] [Google scholar]

Copyright: © 2022 The Authors. This is an open access article under the terms of the Creative Commons Attribution NonCommercial ShareAlike 4.0 (https://creativecommons.org/licenses/by-nc-sa/4.0/). This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.