Research Paper

INTRODUCTION

Global biodiversity is undergoing a precipitous decline, driven by climate change and intensified by anthropogenic pressures (Dirzo and Raven 2003). Among avian species, 1469 are threatened with extinction, including 222 classified as critically endangered (BirdLife International 2018). The White-bellied Heron (WBH, Ardea insignis) stands out as one of the rarest, with fewer than 60 individuals remaining worldwide (Price and Goodman 2015, BirdLife International 2024). Bhutan, a Himalayan nation with a strong conservation ethos, currently harbors the largest known population of 29 individuals recorded in 2025 making it a vital refuge for the species (RSPN 2025).

The distribution of rare species like WBH is shaped by a complex interplay of biotic and abiotic factors (Dendup et al. 2024). Historically, WBH occupied undisturbed wetlands across South and Southeast Asia, including eastern Nepal, the Sikkim Terai, northern Bihar, Bhutan Duars, northern Assam, East Pakistan (now Bangladesh), Arakan, and northern Myanmar (Stanford and Ticehurst 1939, Smythies 1953, Walters 1976, Ali and Ripley 1987, King et al. 2001). Despite this broad historical range, sightings have remained rare and fragmented.

Among Bhutan’s 106 recorded waterbird species (Passang 2018), WBH holds special conservation status under Schedule I of the Forest and Nature Conservation Act, 2023 (RGoB 2023). Its preferred habitat includes mid-elevation riparian zones between 100 and 1500 masl, especially shallow banks with water depths of 30 to 45 cm (RSPN 2011). Seasonal movements are common, with herons shifting between major rivers in winter and smaller tributaries during flooding events (Pradhan 2008: unpublished report, Dorji 2011, RSPN 2011). Breeding individuals may remain within general home ranges year-round and can travel up to 25 kilometers between foraging sites (RSPN 2011, 2015, Price and Goodman 2015).

Nest site selection is influenced by ecological factors such as food availability, predation risk, conspecific presence, nesting material access, and climatic conditions (Collias 1986, Hansell and Hansell 2005, Mainwaring et al. 2014). As a waterbird species, WBH is highly sensitive to habitat degradation, rendering it vulnerable to extinction under accelerating climate change (Maheswaran et al. 2021a). Ensemble modeling reveals a 51.98% loss in WBH habitat which is approximately 12,139 km² across its former range, with the steepest declines in India, Myanmar, and Bangladesh (Maheswaran et al. 2021b). Conversely, Bhutan has experienced a 60% increase in suitable habitat, particularly in northern regions, attributed to intact forest and riverine ecosystems, low human disturbance, and climate-driven shifts in viable zones.

Despite these ecological advantages, Bhutan’s WBH population faces emerging threats. Hydropower development, land-use changes, and climate-induced alterations in river flow are beginning to erode habitat quality (RSPN 2024). The national emphasis on hydropower as a development strategy has led to ecosystem disruptions (Nomura 2025). Additionally, small and fragmented populations remain susceptible to stochastic events such as flash floods and disease outbreaks.

Extensive research has deepened understanding of WBH ecology, including roosting behavior, nest site selection, foraging dynamics, and seasonal movements (Acharja 2019, Khandu et al. 2020a, Khandu et al. 2021). Studies have also explored breeding success, nest predation, and complex social behaviors such as sexual conflict and parental infanticide (Khandu et al. 2020b, Acharja et al. 2021, Khandu 2022). Conservation threats have been assessed through ecological surveys and stakeholder perspectives (Phuntshok et al. 2022, Nima et al. 2025), while acoustic monitoring has provided non-invasive tools for tracking individuals and population trends (Dema et al. 2020). Newly identified critical habitats have expanded known distribution and informed conservation planning (Wangdi et al. 2017).

However, no national-scale assessment has yet isolated the impact of climate change, independent of land-use change on WBH habitat suitability. This gap limits predictive capacity for future conservation planning. The present study addresses this by modeling climate-resilient habitat futures for WBH in Bhutan using coupled model intercomparison project phase 6 (CMIP6) climate scenarios and ensemble modeling. While land-use variables were excluded, the model incorporates bioclimatic predictors, slope, aspect, and proximity to settlements and roads to assess ecological suitability. By quantifying habitat shifts across temporal scales, this study aims to inform adaptive strategies for the long-term survival of this critically endangered species.

METHODOLOGY

Study area

Bhutan is a landlocked Himalayan country in South Asia, covering an area of 38,394 km² and bordered by India to the South, East, and West, and China to the North (27° 31' 53.11" N; 90° 26' 9.07" E). The country experiences four distinct seasons: spring (March–May), summer (June–August), autumn (September–November), and winter (December–February). Situated at the intersection of the Indo-Malayan and Palaearctic realms and with diverse climatic and geographic conditions, Bhutan is part of a global biodiversity hotspot (Myers et al. 2000). It maintains 69.71% forest cover and has placed 52% of its land under formal protection (FMID 2023, WWF 2023).

Bhutan supports rich biological diversity, approximately 200 mammalian species (including tiger, snow leopard, red panda), 766 bird species (including WBH), 5600 plant species, 778 butterfly species, 1940 moth species, and over 130 fish species (Wangchuk et al. 2004, Nepal 2022, Dendup et al. 2023, Dorji et al. 2024, NBC 2025). This diversity is sustained by an extensive protected area (PA) network comprising 11 PAs including the Royal Botanical Park and 9 biological corridors (BC; Fig.1), collectively covering around 52% of the national territory (WWF 2023). The country’s wide elevational gradient, ranging from 97 masl to 7500 masl, supports 11 distinct forest types, from subtropical broadleaf forests to alpine scrub, each offering unique habitat conditions for various species (DoFPS 2022).

Bhutan’s intact forest ecosystems and rugged terrain, dominated by steep mountains, are intersected by four major river basins: Drangmechhu, Punatsangchhu, Amochhu, and Wangchhu (Alam et al. 2017). These river systems provide critical habitat for aquatic fauna and waterbirds, which predominantly occupies these Himalayan rivers. However, long-term monitoring has revealed shifts in WBH distribution, with signs of local extirpation and nesting site abandonment in several historically occupied areas. This trend is particularly pronounced in the Punatsangchhu basin, where upstream sighting frequencies have declined markedly in recent years (RSPN 2025).

Bhutan’s geography and conservation policies continue to offer refuge for WBH and other threatened species, but emerging patterns underscore the need for adaptive management and targeted conservation planning.

Data collection

Species occurrence data collection

For this study, we compiled a total of 615 occurrence records of the WBH from two primary sources. A total of 132 records spanning 2020 to 2024 were obtained from the Spatial Monitoring and Reporting Tool (SMART), an advanced conservation management platform designed to support protected area and wildlife authorities in enhancing patrol monitoring, performance evaluation, and adaptive management strategies (FAO 2025). Field data contributing to SMART are collected by officials from the Department of Forests and Park Services, with each WBH observation precisely geotagged using GPS-enabled devices and timestamped, thereby enabling high-resolution spatial analysis for conservation planning.

The remaining 483 records, covering the period 2001 to 2023, were downloaded from the Global Biodiversity Information Facility (GBIF; https://www.gbif.org/), an internationally supported network and data platform that provides free and open access to biodiversity data from around the world. GBIF enables users to explore information on all forms of life, contributing to global biodiversity research and conservation efforts.

Following a thorough verification process and the removal of duplicate entries (n = 254), we retained 361 unique occurrence points for subsequent analysis.

Environmental variables

Three topographic predictors, elevation, slope, and aspect were initially considered based on their ecological relevance to WBH habitat. These variables were derived from a 12.5-meter resolution Digital Elevation Model (DEM) provided by the Forest Monitoring and Information Division (www.dofps.gov.bt). Elevation can influence WBH occurrence indirectly by shaping temperature, vegetation type, and human accessibility. Slope and aspect were retained, as they may affect river morphology, microclimatic conditions, and riparian vegetation factors that influence prey availability and roosting suitability. This consolidated approach avoids redundancy while clarifying the ecological rationale behind variable selection.

Anthropogenic variables

To account for human influences on WBH distribution, we incorporated two anthropogenic variables: Euclidean distance to settlement (entire settlement of the country - point data) and Euclidean distance to road networks (entire road networks of the country - line data). These spatial datasets were obtained from the National Land Commission Secretariat (https://web.nlcs.gov.bt/). These variables serve as proxies for human accessibility and activity, helping to identify relatively undisturbed riverine stretches preferred by WBH. Their inclusion complements the ecological predictors and enhances the model’s ability to delineate suitable habitat across Bhutan’s landscape.

Bioclimatic variables (projections and downscaling)

Daily gridded temperature and precipitation data at 1 km resolution for the historical period 1996–2014 were obtained from the National Centre for Hydrology and Meteorology (NCHM, https://www.nchm.gov.bt/), and aggregated to monthly resolution. The daily datasets called BhutanClim was prepared using quality-controlled data from 67 weather stations operated by NCHM (Lehner and Formayer 2023). Future climate projections were sourced from ten (10) CMIP6 Global Climate Models (GCMs) under the SSP3 scenario, selected via the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) framework to ensure diversity in model structure, climate sensitivity, and historical performance. The details of the GCMs included is in Table 1.

Monthly climate outputs for two future periods 2041–2060 and 2061–2100 were individually downscaled to 1 km spatial resolution using bilinear interpolation. Bias correction was applied by comparing modeled outputs to BhutanClim observations: temperature was adjusted using mean differences and precipitation using the ratio method. Ensemble means across all ten models were calculated for each climate variable to support ecological impact modeling.

To forecast potential shifts in the distribution of the WBH, climate layers were aligned with the Shared Socioeconomic Pathway 3 (SSP3) scenario, characterized by regional rivalry, limited international cooperation, and significant challenges to mitigation and adaptation (Riahi et al. 2017). This scenario was selected for its conservative representation of a high-emission future. All future projections employed BhutanClim-derived bioclimatic variables at consistent spatial resolution.

Modeling procedures

The modeling of the current and future distribution of the WBH, was based on 361 presence records, supplemented with pseudo-absence points generated at a ratio of 1.5:1. This approach is commonly adopted to improve model discrimination in species distribution modeling (Barbet-Massin et al. 2012). The model building was performed in R v.4.3.3 software (R core team 2024) and further analysis were carried out using ArcMap v.10.8.2 (ESRI 2020)

A total of 24 environmental predictors were selected for habitat suitability modeling, comprising elevation, slope, aspect, Euclidean distances to settlements and road networks, and 19 bioclimatic variables (BIO1-BIO19). The bioclimatic layers were sourced from the BhutanClim dataset at a resolution of approximately 1 kilometer, representing baseline climate conditions from 1996 to 2014.

To harmonize spatial scales and ensure compatibility during modeling, the finer-resolution topographic and anthropogenic data layers were resampled to match the course 1 km resolution of the BhutanClim climate layers. This resampling step was crucial for several reasons: it minimized artificial precision and reduced spatial bias in overlay analyses (Guisan et al. 2007), it ensured scale consistency for algorithms that are sensitive to resolution mismatches (Elith and Leathwick 2009), and it aligned all predictors with the native resolution of climate models, which are typically produced at coarser scales (Peterson et al. 2011).

Before modeling, we applied a multicollinearity filter using a 0.70 correlation threshold, retaining predictors that were ecologically informative yet statistically independent (Dormann et al. 2012). Species distribution modeling was conducted using the sdm package in R (Naimi and Araújo 2016), leveraging an ensemble of Generalized Linear Model (GLM), Generalized Additive Model (GAM), and Random Forest (RF) algorithms. We employed bootstrapping (n = 2) to stabilize model outputs and reduce sampling bias. Model performance was assessed using the True Skill Statistic (TSS), selected over the more common AUC metric due to its robustness for presence-absence datasets (Allouche et al. 2006). Importantly, the ensemble approach combining multiple algorithms minimizes prediction error, improves robustness in spatial probability assessment, and enhances reliability in modeling efforts by combining multiple models into an individual ensemble model (Elith et al. 2006, Franklin 2009, Peterson et al. 2011, Mukherjee et al. 2020, Mohamed et al. 2023, Koç et al. 2024), making them particularly suitable for use in a complex terrain like Bhutan.

Habitat suitability predictions for the WBH were made for two future periods 2041–2060 and 2061–2100 under the SSP3 climate scenario to assess potential habitat conditions in the mid and late 21^st century. Final distribution maps were derived using weighted ensemble outputs combining predictions from GLM, GAM, and RF models.

District-level habitat suitability modeling

Ecological suitability for the WBH across Bhutan was assessed using species distribution modeling and zonal statistics in ArcGIS. In the habitat modeling, a threshold value of 0.5 was used and accordingly habitat rasters were categorized from 0 (unsuitable) to 4 (high suitability) and analyzed across district boundaries to extract metrics such as MEAN, MAJORITY, MAX, and VARIETY. Suitability was evaluated for three time frames namely baseline (1996–2014), 2041–2060, and 2061–2100 under SSP3 climate projections. Temporal comparison of MEAN scores enabled transition analysis, revealing district-specific trends in habitat stability, degradation, or improvement.

Protected area suitability assessment

Habitat suitability within Bhutan’s PA was evaluated using spatial overlays and zonal statistics across three climate scenarios using key metrics as in district-level habitat suitability modeling to assess ecological conditions.

Protection gap and overlay analysis

A protection gap analysis was conducted to identify highly suitable WBH habitats outside Bhutan’s PA. Class 4 habitat zones were isolated using raster reclassification, and overlaps with PA were removed via spatial tools. These unprotected areas were intersected with district boundaries to quantify their extent. A Gap Index was calculated by dividing suitable area outside PA by total district suitability and normalized for comparison.

RESULTS

Multicollinearity screening

Following a Variance Inflation Factor screening to minimize multicollinearity, nine essential predictor variables were selected for modeling. These included topographic features (slope, aspect), anthropogenic factors (Euclidean distance to settlements and roads), and five bioclimatic parameters: BIO3 (isothermality), BIO8 (mean temperature of the wettest quarter), BIO12 (annual precipitation), BIO14 (precipitation of the driest month), and BIO15 (precipitation seasonality). Collectively, these variables capture critical gradients in terrain, human disturbance, and climate across Bhutan, providing a robust framework for predicting the distribution of the WBH (Table 2).

Predictor influence and model performance

To assess the ecological drivers and reliability of the species distribution model, predictor response curves and classification performance metrics were analyzed in tandem. Nine key environmental and anthropogenic variables exhibited threshold-dependent relationships with WBH habitat suitability, with notable avoidance of human-modified landscapes and strong positive associations with slope and annual precipitation. Additionally, all three modeling algorithms GLM, GAM, and RF exhibited varying levels of predictive performance across multiple metrics. While RF and GAM achieved high accuracy (AUC ≥ 0.91), strong correlation (COR ≥ 0.74), and robust skill scores (TSS ≥ 0.74), GLM showed moderate predictive strength (AUC = 0.71, COR = 0.33, TSS = 0.45). Deviance values further highlighted RF as the best-fitting model (0.35), followed by GAM (0.87) and GLM (1.25), reinforcing the advantage of an ensemble approach that integrates complementary strengths (Table 3).

Current habitat suitability (baseline)

The ensemble model predicted 8219 km² of climatically suitable habitat under current conditions. The baseline climate conditions predicted the spatial extent and relative dominance of three habitat suitability classes less suitable, moderately suitable, and highly suitable across the study area. The model identifies 1,985 km² as moderately suitable habitat, with the largest share classified as less suitable (5,305 km²), while highly suitable zones (929 km²) align with core riparian landscapes and mid-elevation regions critical for WBH conservation (Fig. 2).

Future habitat suitability under SSP3 scenarios (2041–2060 and 2061–2100)

Species distribution projections under the SSP3 climate scenario indicate a significant expansion in climatically suitable habitat for the WBH across Bhutan. By mid-century (2041–2060), total suitable habitat increases to approximately 11,980 km² (less suitable 6568; moderately suitable 3376 and highly suitable 2036), marking a net gain of 3761 km² over baseline conditions and the overall change includes 7473 km² of habitat gain and 1464 km² of loss (Fig. 3).

This upward trend continues into the late-century (2061–2100), with suitable habitat reaching 13,784 km² (less suitable 6775 moderately suitable 4202 and highly suitable 2807) a net increase of 5565 km², based on 9955 km² gained and 1128 km² lost (Fig. 4). Despite this growth, the suitability patterns become increasingly fragmented, emphasizing the need for strategies to enhance connectivity and climate-resilient landscape management.

District-level habitat suitability summary: baseline conditions

Under baseline conditions, districts such as Dagana, Tsirang, and Zhemgang demonstrate strong ecological suitability for the WBH, with MAJORITY values of 3 and MEAN scores exceeding 0.75, positioning them as key conservation priorities and potential refugia. Monggar and Pemagatshel exhibit high habitat heterogeneity (VARIETY = 4–5), and although their MEAN scores are moderate (~0.28–0.55), the presence of highly suitable pockets (MAX = 4) suggests strong potential for targeted management. In contrast, districts like Samdrup Jongkhar and parts of Chhukha report MAJORITY values of 1 or 2 and MEAN scores below 0.3, indicating reduced habitat quality likely influenced by anthropogenic pressure or suboptimal terrain. Notably, areas with high VARIETY scores reflect ecotonal transitions and niche diversity, while districts such as Samtse and Trashigang, despite overall low MEAN scores, show high MAX values highlighting the presence of isolated but ecologically valuable habitat fragments deserving strategic attention.

District-level habitat suitability summary for 2041–2060 under SSP3

Under SSP3 climate projections for 2041–2060, Dagana, Pemagatshel, and Monggar maintain high habitat suitability for the WBH, with MEAN scores exceeding 0.6 and dominant suitability classes (MAJORITY = 3 or 4), indicating strong ecological persistence despite projected pressures. Paro, Chhukha, and Bumthang exhibit moderate suitability (MEAN ~0.45–0.50), with habitat classes reaching up to MAX = 3, positioning them as transitional zones where conservation interventions may remain effective. Conversely, Gasa and Haa reflect critically low habitat suitability (MEAN = 0.02) and minimal ecological diversity (VARIETY ≤ 2), likely driven by topographic or climatic limitations (Table 4).

District-level habitat suitability summary for 2061–2100 under SSP3

Under SSP3 projections for 2061–2100, habitat suitability across most districts shows a marked decline compared to baseline and mid-century estimates (Table 5).

Previously resilient zones like Dagana exhibit steep reductions, with MEAN suitability dropping from ~1.16 to ~0.39 and MAJORITY class shifting from 4 to 1, indicating widespread degradation. Moderate retention zones such as Pemagatshel, Monggar, and Tsirang sustain MEAN values between ~0.35 and 0.43, suggesting the presence of fragmented but viable habitat patches. In contrast, critically low suitability persists in districts like Gasa, Haa, Wangduephodrang, and Bumthang, where MEAN scores fall below 0.10 and MAJORITY suitability is absent, likely driven by climatic shifts, topographic constraints, or increasing anthropogenic pressures (Table 6).

Transition insights

Between 2041 and 2100, modeled habitat suitability trends reveal three distinct trajectories across districts. Dagana and Pemagatshel initially show gains in ecological quality during mid-century projections, but subsequently decline sharply, likely reflecting early conservation benefits offset by intensified SSP3 stressors. In contrast, districts such as Gasa, Haa, and Trashigang exhibit consistent deterioration throughout both periods, ending with minimal suitability and limited conservation potential. Notably, Monggar demonstrates gradual improvement over time, resulting in a net gain in habitat suitability by 2100. This trajectory positions Monggar as a potential future refugium and target for climate-resilient conservation planning (Table 7).

Current habitat suitability within PA (species distribution modeling baseline scenario)

Under baseline conditions, Royal Manas National Park (RMNP) and Phibsoo Wildlife Sanctuary (PWS) stand out as high-performing PA, exhibiting MEAN suitability scores between 0.9 and 1.2 and dominant habitat classes (MAJORITY = 3 or 4). These sites serve as critical ecological strongholds for biodiversity conservation. In contrast, Jomotsangkha Wildlife Sanctuary demonstrates moderate habitat suitability (MEAN = 0.4-0.7) alongside high internal diversity (VARIETY = 4–5), suggesting its potential to support species persistence if adequately buffered. Meanwhile, high-altitude or rugged terrain PA such as Wangchuck Centennial Park and Bumdeling Wildlife Sanctuary (BWS) show low suitability scores (MEAN < 0.3; MAJORITY = 1–2), likely reflecting climatic constraints or topographic barriers to habitat quality.

Projected habitat suitability within PAs (2041–2060, SSP3 scenario)

Under projected SSP3 conditions for 2061–2100, distinct habitat suitability patterns emerged across Bhutan’s PA. Southern corridors BC4, BC3, and BC2 demonstrated strong ecological resilience, with MEAN scores ranging from 0.78 to 0.95 and dominant suitability classes (MAJORITY = 3). Their internal habitat diversity (VARIETY = 4) underscores their importance as long-term refugia for the WBH.

Moderate zones such as BC5, BC7, and Jigme Singye Wangchuck National Park (JSWNP) displayed intermediate MEAN scores (~0.30–0.45) and mixed suitability classes, indicating fragmented but potentially restorable landscapes. These areas may benefit from strategic buffering or connectivity enhancement to reinforce ecological viability.

In contrast, BC6, Jigme Khesar Strict Nature Reserve (JKSNR), and BWS reported severe ecological decline, with MEAN values ≤ 0.01, low maximum suitability (MAX ≤ 1), and MAJORITY classes of 0. These zones are likely unsuitable for future WBH persistence and may require alternate approaches, such as elevational rewilding or long-term restoration planning.

Projected habitat suitability within PAs (2061–2100, SSP3 scenario)

Under the SSP3 scenario for 2061–2100, PA zones exhibited clear contrasts in habitat quality for the WBH. Southern corridors namely BC2, BC3, and BC4 retained strong ecological integrity, with MEAN scores exceeding 1.07 and high habitat suitability classes persisting. These areas remain critical strongholds for biodiversity under intensifying climate pressures.

Moderate-performing zones such as BC5, BC7, and JSWNP exhibited viable habitat conditions, but showed notable internal variability in suitability scores, as reflected by standard deviation values exceeding 1.0. This high STD indicates spatial fragmentation within these zones i.e., a mix of highly suitable and unsuitable habitat patches. Such heterogeneity may affect habitat continuity and could warrant targeted micro-zonation and internal restoration planning to enhance ecological connectivity.

Conversely, BC6, BWS, JKSNR, and Jigme Dorji National Park (JDNP) displayed severe habitat decline, with MEAN values as low as 0.001-0.11 and suitability capped at Class 1. Their poor resilience under SSP3 suggests a need for alternative conservation strategies, such as elevational rewilding or assisted migration planning.

Consolidated habitat suitability findings in districts and PAs (baseline → 2041-2060 → 2061-2100 SSP3)

Eastern districts, particularly Monggar and Pemagatshel, consistently retain areas of suitable WBH habitat that lie outside Bhutan’s PA network across all modeled timeframes. These districts are strong candidates for recognition as Other Effective Area-Based Conservation Measures (OECMs).

Mid-elevation regions adjoining the BC3 and BC4 corridors consistently demonstrate ecological resilience across both baseline and future SSP3 scenarios, underscoring their importance in climate-adaptive conservation strategies for the WBH. In contrast, northern zones such as JKSNR have never been identified as suitable habitat, while BWS shows low suitability for the baseline, but persistent unsuitability for future timeframes. These findings highlight the limited ecological viability of high-altitude areas and reinforce the need to redirect conservation planning toward regions with demonstrated habitat potential and long-term resilience. (Table 8).

Protection gap analysis

Distribution of suitable WBH habitat outside PA (2041–2060, 2061–2100)

Across most districts, only a very small proportion of suitable WBH habitat falls outside PA, typically ranging from 0.0% to 1.5%, with slight increases in fragmentation projected by 2061–2100. Districts like Lhuentse, Monggar, Paro, and Wangduephodrang consistently show moderate gaps (~0.7–0.9%) throughout both future periods, indicating the presence of minor buffer zones or transitional fragments that may require further protection. Trongsa exhibits a more pronounced rise in unprotected habitat—from 0.3% to 1.5%—signaling potential fragmentation pressure in central Bhutan that could benefit from ecological corridor planning. Meanwhile, districts such as Gasa, Haa, Pemagatshel, Samtse, Trashigang, and Yangtse report no suitable WBH habitat outside PA (Table 9).

DISCUSSION

The results of this study provide the most comprehensive climate-informed distribution modeling for the WBH in Bhutan to date, with significant implications for conservation planning in the eastern Himalayas. Ensemble projections across SSP3 scenarios present compelling evidence that, while Bhutan currently serves as a critical refugium, future climate trajectories are likely to reshape habitat distributions, fragmenting current strongholds and introducing novel, often unprotected, high-suitability zones. These findings underscore the urgent need for adaptive, anticipatory conservation strategies rooted in spatially explicit data and long-term ecological forecasting.

Climate-driven range expansion and retraction

Contrary to negative trends observed in other range countries such as India, Myanmar, and Bangladesh (Maheswaran et al. 2021b), Bhutan demonstrates a significant net habitat gain for WBH across all projected timeframes. Climatically suitable areas expand from 8219 km² under baseline conditions to 11,980 km² by 2041–2060, and further to 13,784 km² by 2100. This trajectory indicates that Bhutan’s mid-elevation ecosystems, particularly riparian zones buffered by intact forest cover, may become increasingly vital for global WBH conservation echoing broader range-shift patterns in Himalayan taxa under thermal stress (Telwala et al. 2013).

However, expansion is neither linear nor universally beneficial. Notably, high-suitability areas begin to emerge in previously low/unreported WBH presences districts such as Monggar, Pemagatshel, and Lhuentse, many of which lie outside the existing PA network. These new zones mirror range shifts observed in other montane bird species, where climate change drives habitat redistribution beyond current conservation boundaries (Hole et al. 2009). This emerging mismatch between climate-adjusted habitat distributions and formal conservation coverage demands recalibration of Bhutan’s spatial planning.

Fragmentation and ecological instability

Despite aggregate habitat gains, model outputs suggest an increase in spatial fragmentation. Districts such as Dagana and Tsirang, which demonstrate strong suitability in mid-century projections, undergo steep declines by 2100. Northern zones, including Gasa and Haa, remain persistently unsuitable across all scenarios, likely constrained by topography, altitude, and climatic variability.

This fragmentation reduces carrying capacity and amplifies risks for isolated WBH populations already vulnerable due to low reproductive rates and susceptibility to stochastic events like flooding or disease outbreaks (Price and Goodman 2015). The emergence of high-suitability patches surrounded by suboptimal landscapes reflects increasing ecological isolation, a trend consistent with climate-induced fragmentation documented among river-dependent species (Domisch et al. 2013). Moreover, PA such as JSWNP, BC5, and BC7 show moderate MEAN suitability but high internal variability (STD > 1), emphasizing the need for micro-zonation and habitat restoration to preserve internal connectivity.

Protected area efficacy and coverage gaps

One of the most striking results is the mismatch between projected WBH habitat and Bhutan’s existing conservation network. While core PA like RMNP and PWS continue to act as ecological anchors, emerging high-suitability zones in districts such as Monggar, Trongsa, and Dagana fall outside PA boundaries.

The gap analysis reveals that districts like Tsirang, Paro, and Wangduephodrang may have over 80% of their future WBH habitat outside formal conservation areas. These spatial gaps highlight the need to move beyond a static PA model toward a more dynamic, climate-responsive network.

The integration of Other Effective Area-Based Conservation Measures (OECMs) offers a promising strategy for biodiversity conservation in Bhutan. Recognized under the Convention on Biological Diversity’s post-2020 global framework, OECMs encompass community-managed forests, sacred sites, riparian buffer zones, and seasonal no-development zones that deliver long-term biodiversity outcomes (IUCN-WCPA 2019). Districts like Monggar and Pemagatshel, consistently retaining suitable habitat outside the PA system, are strong candidates for OECM designation. This approach allows Bhutan to formalize informal stewardship practices while enhancing ecological coverage.

Conservation planning in a fragmented future

The study’s findings offer actionable guidance for Bhutan’s conservation planning. Protected corridors such as BC2, BC3, and BC4 demonstrate high habitat retention across all time frames, positioning them as potential climate-resilient refugia and critical dispersal routes. Maintaining and enhancing landscape connectivity between these corridors especially across elevational gradients should become central to WBH conservation policy.

Conversely, JKSNR was never identified as suitable habitat for WBH in our model; it remains ecologically irrelevant to the species. Although BWS recorded a single vagrant sighting, it was not modeled as suitable habitat. Thus, references to suitability loss in these areas do not indicate future decline but rather reflect their persistent unsuitability. In contrast, both current and future suitable habitats are concentrated along BC8, particularly in the Kurichu basin; however, these areas are not currently occupied by WBH. This distinction between modeled suitability and actual occupancy is essential for guiding conservation priorities and avoiding misallocation of resources. Moderate-performing zones like JSWNP and BC7, though spatially fragmented, can benefit from targeted internal management interventions e.g., riparian reforestation, or seasonal visitor zoning.

Importantly, districts containing small but high-quality habitat pockets evidenced by high MAX scores amidst declining MEAN values may serve as ecological stepping stones or satellite habitats. These areas warrant protection through either PA expansion or OECM recognition to ensure spatial continuity.

Wildlife tourism as a conservation ally

In addition to conventional protection mechanisms, wildlife tourism offers a promising yet underutilized tool for WBH conservation in Bhutan. WBH watching activities, designed for non-breeding seasons and developed in coordination with communities and park authorities, could support low-impact ecotourism while fostering public awareness and local stewardship.

Bhutan’s success in nature-based tourism e.g., snow leopard treks, red panda watch in JDNP provides a viable model (Dendup et al. 2021). By establishing observation platforms, interpretive trails, and seasonal guidelines, future WBH tourism initiatives could serve both conservation and livelihoods, especially in emerging refugia districts like Dagana, Monggar, and Pemagatshel. However, careful zoning and impact monitoring will be crucial to prevent disturbance or over-commercialization of sensitive sites.

Policy and research implications

The projected range shifts under SSP3 carry deep policy relevance. Bhutan’s constitutional commitment to maintaining ≥60% forest cover (RGoB 2008) provides a stable foundation for integrating species distribution modeling outputs into national land-use plans, Environmental Impact Assessments, and energy infrastructure development. Given the documented impacts of hydropower on WBH habitat e.g., Punatsangchhu basin degradation strategies to reconcile energy goals with biodiversity priorities are vital (Phuntshok et al. 2022).

On the research front, this study fills a longstanding gap by offering the first nationwide, ensemble-modeled WBH habitat forecast using robust predictors and multitemporal SSP3 projections. Future efforts should focus on validating these predictions through field surveys in eastern and central Bhutan, where model-identified zones have been historically under-sampled. Moreover, incorporating variables like river hydrology, prey availability, and seasonal habitat dynamics can further refine species distribution models for climate-vulnerable taxa. These applications must be tempered by an understanding of the species’ biological constraints and the dynamic nature of anthropogenic pressures, as discussed in the following section.

Limitations

Forecasting habitat suitability for the WBH, a species with fewer than 60 individuals and limited dispersal capacity, presents inherent uncertainties. The species’ small population size, restricted range, and sensitivity to disturbance mean that even minor landscape changes can have disproportionately large ecological impacts. While our models offer valuable insights into potential climate-driven habitat shifts through 2100, they do not account for dynamic anthropogenic pressures such as hydropower development, road expansion, and land-use change, which are expected to intensify in Bhutan over the coming decades. These factors could significantly alter habitat quality and connectivity, regardless of climatic suitability.

Our modeling framework did not explicitly incorporate spatial data on existing or planned hydropower dam sites due to limited access to publicly available, georeferenced datasets. Although Euclidean distance to settlement points was used as a proxy for anthropogenic pressure, this approach may not fully capture the localized impacts of infrastructure such as PHPA I and II, or proposed projects including Sunkosh, Dorjilung, and Chamkharchu. As a result, areas identified as climate-resilient refugia, including parts of Pemagatshel, Dagana, and lower Punatsangchu, may be subject to anthropogenic pressures that were not modeled. Future assessments that integrate infrastructure footprints more directly will be essential to improve predictive accuracy and evaluate long-term habitat viability for WBH.

Additionally, the modeling framework operates at a landscape scale and may overestimate suitable habitat by not fully capturing fine-scale ecological and behavioral filters. For example, predicted suitability in regions such as Paro and northern Bhutan includes areas where WBH has been recorded only rarely. This limited presence is likely due to micro-habitat constraints such as riverbank structure, prey availability, and disturbance intensity. Although northern zones like Gasa and Haa remain persistently unsuitable across all scenarios, this unsuitability is likely driven by topographic and climatic variability. WBH presence is typically restricted to specific stretches of river within broader zones that appear suitable at the landscape scale. We acknowledge this limitation and propose that future research incorporate micro-habitat assessments within confirmed WBH territories to refine model accuracy and enhance ecological realism.

Therefore, our results should be interpreted as indicative rather than prescriptive. They are intended to guide adaptive, precautionary conservation planning that integrates both climate projections and localized habitat conditions.

CONCLUSION

Bhutan’s projected expansion of climatically suitable WBH habitat represents a rare conservation opportunity amidst widespread species decline. However, spatial fragmentation, protection gaps, and ecological volatility underscore the need for multifaceted conservation strategies. By integrating climate-resilient corridor planning, OECMs recognition, regulated wildlife tourism, and fine-scale ecological management, Bhutan can set a global example of adaptive conservation grounded in cultural values and scientific foresight.

Ethic statement

The data used in this study (especially from SMART) were used under the permission of Department of Forests and Park Services.

Declaration of competing interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Funding

No separate funding was received for this work except for the training program attended by the authors (excluding Phub Dhendup and Pema Wangda).

LITERATURE CITED

Acharja, I. P. 2019. Evaluation of nest habitat, site preferences and architecture of the critically endangered White-bellied Heron (Ardea insignis) in Bhutan. Bird Conservation International 30(4):599-617. https://doi.org/10.1017/S095927091900042X

Acharja, I. P., S. Tshering, T. Lhendup, and T. Phuntsho. 2021. First observation of sexual conflict and parental infanticide in the critically endangered White-bellied Heron (Ardea insignis). Birding Asia 35:38-42.

Alam, F., Q. Alam, S. Reza, S. Khurshid-Ul-Alam, K. Saleque, and H. Chowdhury. 2017. Sourcing green power in Bhutan: a review. Energy Procedia 110:586-591. https://doi.org/10.1016/j.egypro.2017.03.189

Ali, S., and S. D. Ripley. 1987. Compact handbook of the birds of India and Pakistan. Oxford University Press, Delhi, India.

Allouche, O., A. Tsoar, and R. Kadmon. 2006. Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic. Journal of Applied Ecology 43(6):1223-1232. https://doi.org/10.1111/j.1365-2664.2006.01214.x

Barbet-Massin, M., F. Jiguet, C. H. Albert, and W. Thuiller. 2012. Selecting pseudo-absences for species distribution models: how, where and how many? Methods in Ecology and Evolution 3(2):327-338. https://doi.org/10.1111/j.2041-210X.2011.00172.x

BirdLife International. 2018. State of the world’s birds: taking the pulse of the planet. BirdLife International, Cambridge, UK.

BirdLife International. 2024. Species factsheet: Ardea insignis (White-bellied Heron). BirdLife International, Cambridge, UK. IUCN Red List 2024. https://www.iucnredlist.org/fr/search/map?taxonomies=22672843&searchType=species

Collias, N. E. 1986. Engineering aspects of nest building by birds. Endeavour 10(1):9-16. https://doi.org/10.1016/0160-9327(86)90044-X

Dema, T., M. Towsey, S. Sherub, J. Sonam, K. Kinley, A. Truskinger, M. Brereton, and P. Roe. 2020. Acoustic detection and acoustic habitat characterization of the critically endangered White-bellied Heron (Ardea insignis) in Bhutan. Freshwater Biology 65(1):153-164. https://doi.org/10.1111/fwb.13217

Dendup, P., R. Dorji, Y. Jamtsho, W. Wangchuk, B. Tshering, and R. Dorji. 2021. Conservation Management Plan of Jigme Dorji National Park for the period July 2021-June 2031. Biodiversity Conservation in Pursuit of Gross National Park. Department of Forests and Park Services, Ministry of Agriculture and Forests. Damji, Gasa, Bhutan.

Dendup, P., P. Kuenzang, T. Dorji, and W. Wangchuk. 2024. Biophysical characteristics, distribution and conservation challenges of critically endangered Bhutan poppy (Papaver bhutanicum) in Jigme Dorji National Park, Bhutan. Integrative Conservation 3(1):76-86. https://doi.org/10.1002/inc3.40

Dendup, P., C. Lham, W. Wangchuk, and Y. Jamtsho. 2023. Tiger abundance and ecology in Jigme Dorji National Park, Bhutan. Global Ecology and Conservation 42:e02378. https://doi.org/10.1016/j.gecco.2023.e02378

Department of Forests and Park Services (DoFPS). 2022. Forest Types of Bhutan. Maps and Statistics. Forest Resources Management Division, Department of Forests and Park Services, Ministry of Agriculture and Forests, Thimphu, Bhutan.

Dirzo, R., and P. H. Raven. 2003. Global state of biodiversity and loss. Annual Review of Environment and Resources 28:137-167. https://doi.org/10.1146/annurev.energy.28.050302.105532

Domisch, S., M. B. Araújo, and D. Eme. 2013. Predicting freshwater fish distributions under climate change: consequences for climate impact assessments. Diversity and Distributions 19(6):837-848.

Dorji, J. 2011. Protecting White-Bellied Heron habitat: extent of anthropogenic threats and people’s attitude towards their conservation in Bhutan. MSc dissertation, College of Natural Resources, Lobesa, Bhutan.

Dorji, R., Letro, S. Yangden, P. Dendup, T. Dhendup, and Y. Lhamo. 2024. Rare and unusual snow leopard encounters in the broadleaf forest of the Bhutanese Himalayas. Snow Leopard Reports 3:13-20. https://doi.org/10.56510/slr.v3.24118

Dormann, C. F., J. Elith, S. Bacher, C. Buchmann, G. Carl, G. Carré, J. R. G. Marquéz, B. Gruber, B. Lafourcade, P. J. Leitão, T. Münkemüller, C. McClean, P. E. Osborne, B. Reineking, B. Schröder, A. K. Skidmore, D. Zurell, and S. Lautenbach. 2012. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36(1):27-46. https://doi.org/10.1111/j.1600-0587.2012.07348.x

Elith, J., and J. R. Leathwick. 2009. Species distribution models: ecological explanation and prediction across space and time. Annual Review of Ecology, Evolution, and Systematics 40:677-697. https://doi.org/10.1146/annurev.ecolsys.110308.120159

Elith, J., C. H. Graham, R. P. Anderson, M. Dudík, S. Ferrier, A. Guisan, R. J. Hijmans, F. Huettmann, J. R. Leathwick, A. Lehmann, J. Li, L. G. Lohmann, B. A. Loiselle, G. Manion, C. Moritz, M. Nakamura, Y. Nakazawa, J. M. M. Overton, A. T. Peterson, and E. E. Zimmermann. 2006. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29(2):129-151. https://doi.org/10.1111/j.2006.0906-7590.04596.x

Environmental Systems Research Institute (ESRI). 2020. ArcGIS Desktop: Release 10.8.2. Redlands, CA, USA.

Food and Agriculture Organization (FAO). 2025. Spatial monitoring and reporting tool (SMART). https://www.fao.org/sustainable-forest-management/toolbox/tools/tool-detail/en/c/1255907/

Forest Monitoring and Information Division (FMID). 2023. National Forest Inventory Volume I: State of Forest Report. Department of Forests and Park Services, Thimphu, Bhutan.

Franklin, J. 2009. Mapping species distributions: spatial inference and prediction. Cambridge University Press, Cambridge, UK. https://doi.org/10.1017/CBO9780511810602

Guisan, A., W. Thuiller, and N. E. Zimmermann. 2007. Habitat modelling using ensemble forecasts. Ecological Modelling 199(1):117-125.

Hansell, M., and M. H. Hansell. 2005. Animal architecture. Oxford University Press, Oxford, UK. https://doi.org/10.1093/acprof:oso/9780198507529.001.0001

Hole, D. G., S. G. Willis, D. J. Pain, L. D. C. Fishpool, S. H. M. Butchart, Y. C. Collingham, and B. Huntley. 2009. Projected impacts of climate change on a continent-wide protected area network. Ecology Letters 12(5):420-431. https://doi.org/10.1111/j.1461-0248.2009.01297.x

IUCN-WCPA. 2019. Guidelines for recognizing other effective area-based conservation measures. IUCN, Gland, Switzerland.

King, B. F., H. Nickel, W. Suter, J. C. Eames, C. Robson, C. Zöckler, and A. Bräunlich. 2001. Birds recorded during two expeditions to North Myanmar (Burma). Forktail 17:1-15.

Khandu, P., G. A. Gale, K. Kinley, T. Tandin, S. Shimano, and S. Bumrungsri. 2020a. Daily roosting behaviour of the critically endangered White-bellied Heron Ardea insignis as a function of day length. Biological Rhythm Research 53(5):812-822. https://doi.org/10.1080/09291016.2020.1814525

Khandu, P., G. A. Gale, R. Pradhen, I. P. Acharya, and S. Bumrungsri. 2020b. First record of successful breeding of the critically endangered White-bellied Heron (Ardea insignis) in broadleaved trees. The Journal of Animal & Plant Sciences 30(2):502-507. https://doi.org/10.36899/JAPS.2020.2.0056

Khandu, D., S. Dorji, and Y. Phuntsho. 2021. Anthropogenic threats to riverine birds in Bhutan. Avian Conservation Bulletin 12(4):203-215.

Khandu, P. 2022. Nest predation of the critically endangered White-bellied Heron (Ardea insignis) by Masked Palm Civet (Paguma larvata) in Burichhu, Bhutan. Journal of Heron Biology and Conservation 7:4.

Koç, D. E., B. Ustaoğlu, and D. Biltekin. 2024. Effect of climate change on the habitat suitability of the relict species Zelkova carpinifolia using ensemble species distribution modelling. Scientific Reports 14:5618. https://doi.org/10.1038/s41598-024-78733-4

Lehner, F., and H. Formayer. 2023. Insights on the climate of Bhutan from a new daily 1 km gridded data set for temperature and precipitation. International Journal of Climatology 43(11):4927-4943.

Maheswaran, B., A. Umesh, and M. Ramkumar. 2021a. Status and distribution of the White-bellied Heron in the Himalayan region. BirdLife Asia Report Series 22:19-36.

Maheswaran, G., L. K. Sharma, H. S. Mondal, and T. Mukherjee. 2021b. White-bellied heron a species on the verge of extinction: ensemble model reveals loss of habitats and resultant prolonged isolation driving the species to extinction. Ecological Informatics 64:101383. https://doi.org/10.1016/j.ecoinf.2021.101383

Mainwaring, M. C., I. R. Hartley, M. M. Lambrechts, and D. C. Deeming. 2014. The design and function of birds’ nests. Ecology and Evolution 4(20):3909-3928. https://doi.org/10.1002/ece3.1054

Myers, N., R. A. Mittermeier, C. G. Mittermeier, G. A. B. da Fonseca, and J. Kent. 2000. Biodiversity hotspots for conservation priorities. Nature 403(6772):853-858. https://doi.org/10.1038/35002501

Mohamed Nisin, K. M. N., K. R. Sreenath, and M. P. Sreeram. 2023. Change in habitat suitability of the invasive Snowflake coral (Carijoariisei) during climate change: an ensemble modelling approach. Ecological Informatics 76:102145. https://doi.org/10.1016/j.ecoinf.2023.102145

Mukherjee, T., L. K. Sharma, G. K. Saha, and K. Chandra. 2020. Past, present and future: combining habitat suitability and future landcover simulation for long-term conservation management of Indian rhino. Scientific Reports 10:606. https://doi.org/10.1038/s41598-020-57547-0

Naimi, B., and M. B. Araújo. 2016. sdm: a reproducible and extensible R platform for species distribution modelling. Ecography 39(4):368-375. https://doi.org/10.1111/ecog.01881

National Biodiversity Centre (NBC). 2025. Species diversity of Bhutan. https://nbc.gov.bt/species-diversity/

Nepal, T. K. 2022. An overview of biodiversity in Bhutan. Asian Journal of Research in Agriculture and Forestry 8:7-19. https://doi.org/10.9734/ajraf/2022/v8i130144

Nima, P., T. Dorji, M. S. Rana, and T. Dorji. 2025. Knowledge, attitudes, perceived threats, and conservation challenges of the critically endangered White-bellied Heron (Ardea insignis) in Bhutan. Global Ecology and Conservation 58:e03484. https://doi.org/10.1016/j.gecco.2025.e03484

Nomura, K. 2025. Bhutan’s hydropower-led economic growth. Pages 1-13 in Hydropower-led economic growth in Bhutan. Contemporary South Asian Studies. Springer Nature, Singapore. https://doi.org/10.1007/978-981-96-5150-4_1

Passang. 2018. Diversity of waterbirds along Bindu River of Tashi Chhoeling (Samtse), Bhutan. Journal of the Bhutan Ecological Society 3:45-55.

Peterson, A. T., J. Soberón, R. G. Pearson, R. P. Anderson, E. Martínez-Meyer, M. Nakamura, and M. B. Araújo. 2011. Ecological niches and geographic distributions. Princeton University Press, Princeton, New Jersey, USA. https://doi.org/10.23943/princeton/9780691136868.003.0003

Phuntshok, P., K. Sherub, and N. Dorji. 2022. Conservation threats to White-bellied Heron (Ardea insignis) in Zhemgang District, Bhutan. Asian Journal of Research in Zoology 5(1):14-20. https://doi.org/10.9734/ajriz/2022/v5i130127

Price, M. R. S., and G. L. Goodman. 2015. White-bellied Heron conservation strategy 2015. Royal Society for Protection of Nature. Thimphu, Bhutan.

R Core Team. 2024. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://doi.org/10.32614/R.manuals

Riahi, K., D. P. Van Vuuren, E. Kriegler, J. Edmonds, B. C. O’Neill, S. Fujimori, N. Bauer, K. Calvin, R. Dellink, O. Fricko, W. Lutz, A. Popp, J. C. Cuaresma, S. Kc, M. Leimbach, L. Jiang, T. Kram, S. Rao, J. Emmerling, and M. Tavoni. 2017. The shared socioeconomic pathways and their energy, land use, and greenhouse gas emissions implications: an overview. Global Environmental Change 42:153-168. https://doi.org/10.1016/j.gloenvcha.2016.05.009

Royal Government of Bhutan (RGoB). 2008. Constitution of the Kingdom of Bhutan. Government Press. Thimphu, Bhutan.

Royal Government of Bhutan (RGoB). 2023. Forest and Nature Conservation Act of Bhutan, 2023. Thimphu, Bhutan.

Royal Society for Protection of Nature (RSPN). 2011. The critically endangered White bellied Heron. R. Pradhan, and P. Frederick, editors. Royal Society for Protection of Nature. Thimphu, Bhutan.

Royal Society for Protection of Nature (RSPN). 2015. Strengthening conservation of critically endangered White-bellied Heron. Rangzhin 8(4):1-2.

Royal Society for Protection of Nature (RSPN). 2024. White-bellied heron (Ardea insignis) annual population survey 2024. Royal Society for Protection of Nature. Thimphu, Bhutan.

Royal Society for Protection of Nature (RSPN). 2025. Annual White-bellied Heron population survey report. Royal Society for Protection of Nature. Thimphu, Bhutan.

Smythies, B. E. 1953. Birds of Burma. 2nd edition. Oliver & Boyd, Edinburgh, UK.

Stanford, J. K., and C. B. Ticehurst. 1939. On the birds of northern Burma. Ibis Series 14(2):65-638. https://doi.org/10.1111/j.1474-919X.1939.tb03970.x

Telwala, Y., B. W. Brook, K. Manish, and M. K. Pandit. 2013. Climate-induced elevational range shifts and increase in plant species richness in a Himalayan biodiversity epicentre. PLoS ONE 8(6):e68179. https://doi.org/10.1371/journal.pone.0057103

Walters, M. 1976. Some observations on the eggs of the Great White-bellied Heron, Ardea insignis. Journal of the Bombay Natural History Society 73(1):213-214.

Wangchuk, T., P. Thinley, K. Tshering, C. Tshering, D. Yonton, and B. Pema. 2004. Mammals of Bhutan. Department of Forests and Park Services, Thimphu, Bhutan.

Wangdi, K., T. Dhendup, and T. Tshering. 2017. A new critical habitat for conservation of the White-bellied Heron Ardea insignis Hume, 1878 (Aves: Ardeidae) from Bhutan. Journal of Threatened Taxa 9(2):9862 -9863. https://doi.org/10.11609/jott.2892.9.2.9862-9863

World Wildlife Fund (WWF). 2023. Biological Corridor - 9: a green connecting Sakteng Wildlife Sanctuary and Bumdeling Wildlife Sanctuary, Trashiyangtse, Bhutan. https://www.wwfbhutan.org.bt/?381915/Biological-corridor-9#:~:text=This%20increases%20Bhutan's%20protected%20area,Bhutan's%20conservation%20journey%20till%20date