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

Research Article - Journal of Environmental and Occupational Health (2022)

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

Sangay Dorji1*, Lakey Lakey2, Tshering Wangchen2 and Steve Adkins1
1Department of Agriculture and Food Sciences, The University of Queensland, Brisbane, Australia
2Department of Agriculture, Ministry of Agriculture and Forests, Thimphu, Bhutan
*Corresponding Author:

Sangay Dorji, Department of Agriculture and Food Sciences, The University of Queensland, Brisbane, Australia, Email:

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


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 km2) 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 km2 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.


Climate change; Biomod; Parthenium weed; Species distribution modelling


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 km2, 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 ( [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 ( [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 km2) 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].


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).

Table 1: Performance score of individual model based on the Receiver Operating Characteristics (ROC) and the Area Under the Curve (AUC) of and the True Skills Statistics (TSS). AUC and TSS presented is the average of five evaluation runs and three pseudo-absences run of each modelling algorithm.


















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.

Table 2: Importance of predictor variables used in distribution modelling of P. hysterophorus in Bhutan. Variable importance is presented as the mean over five evaluation runs and three pseudo-absences datasets of each modelling algorithm, and as the mean of means amongst them. The models: Generalized Additive Model (GAM), Generalized Boosted Model (GBM), Generalized Linear Model (GLM), Random Forest (RF), Maximum Entropy (MAXENT.Phillips). Variables of greatest influence are highlighted in bold.

Model Aspect *_Bio15 *_Bio18 *_Bio19 *_Bio2 *_Bio3 *_Bio6 *_Bio8 Slope


0.1353 0.5723 0.7856 0.2431 0.2253 0.2691 0.3478 0.8304 0.0675


0.0285 0.0993 0.0142 0.2921 0.0320 0.0039 0.0069 0.6434 0.0027


0.1277 0.4938 0.3348 0.1929 0.1383 0.1040 0.3701 0.9205 0.0199


0.0149 0.0379 0.0240 0.1437 0.0352 0.0097 0.0701 0.2038 0.0039


0.1937 0.2331 0.4140 0.1498 0.1659 0.0543 0.6685 0.0685 0.1455

Mean of means

0.1000 0.2873 0.3145 0.2043 0.1193 0.0882 0.2927 0.5333 0.0479

Note: *Bhutan


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 km2 (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 km2 (12.98%), 129.31 km2 (11.76%) and 107.01 km2 (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

Table 3: Current and future predicted habitat suitability areas in Bhutan

Period Suitability (Km2) Area change (Km2) Percent (%) of current predicted distribution Percent (%) of Bhutan’s area Status
Current 1099.10     2.83  
RCP 4.5 2050 5419.70 4320.60 393.10 13.95 Gain
RCP 8.5 2050 3364.18 2265.08 206.09 8.66 Gain
RCP 4.5 2070 2870.37 1771.27 161.16 7.39 Gain
RCP 8.5 2070 5054.07 3954.97 359.84 13.01 Gain

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 km2 (13.95% of Bhutan’s area), followed by RCP8.5 2070 with about 3,954.97 km2 (13.01% of Bhutan’s area). The lowest is predicted under RCP4.5 2070 with about 1,771.27 km2 (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 km2 followed by Dagana with about 692.23 km2. Under RCP8.5 2050 scenario, Dagana had the highest suitability with about 462.60 km2 followed by Wangduephodrang with about 429 km2. Under RCP4.5 2070 scenario, Dagana had the highest suitability with about 441.42 km2 followed by Wangduephodrang with about 434.74 km2. Under RCP8.5 2070 scenario, the highest suitability was predicted in Dagana with about 681.08 km2 followed by Wangduephodrang with about 560.69 km2.

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 km2 while the lowest was determined under RCP 4.5 2070 scenario with about 1,959.65 km2 (Figure 4 and Table 4). In terms of loss, the highest loss in suitability was determined under RCP4.5 2070 with about 178.35 km2 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 km2 while the lowest was predicted under RCP4.5 2070 with about 957.53 km2 (Figure 4 and Table 4).

Table 4: Change in the distribution of P. hysterophorus between current and future periods.

Period Loss Stable Gain
RCP 4.5 2050 5.57 1130.31 4391.94
RCP 8.5 2050 27.87 1108.02 2334.19
RCP 4.5 2070 178.35 957.53 1959.65
RCP 8.5 2070 82.49 1053.40 4051.95

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.


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.


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.


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.


Copyright: © 2022 The Authors. This is an open access article under the terms of the Creative Commons Attribution NonCommercial ShareAlike 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.