The great hornbill (Buceros bicornis) and rufous-necked hornbill (Aceros nepalensis) are listed as vulnerable under the International Union for Conservation of Nature Red List of threatened species due to the rapid decline of their population in the world. This research focuses on analyzing the habitat suitability of these two important bird species across Bhutan. A total of 51 presence locations were recorded from the field survey. The models were simulated using three topographic variables and 19 bioclimatic variables. The MaxEnt modelling technique was used for delineating the distribution potential habitat suitability map. The habitat suitability analysis for great hornbill and rufous-necked hornbill shows that 2% and 3% of Bhutan’s total geographical area were highly suitable, respectively. The approach of this study will be beneficial in identifying suitable areas and aid decision-makers in management and conservation of these vulnerable bird species.

1. Introduction

Species distribution modelling (SDM) is a widely used methodology that has been used to predict the distribution of species on any geographical range [1]. These are also known as bioclimatic models envelope models and ecological niche models [2]. The Maximum Entropy (MaxEnt) model predicts the potential habitat considering entropy of different variables associated with present location of species [3]. SDM is one of the most active areas of global ecology, and several papers have been published around the world to date [4]. These are models that use species environment relationships to explain and predict distributions of species [5]. They rely on statistical correlations between existing species distributions and environmental variable [6].

Hornbills are recognized as one of the most important bird species in the tropical and subtropical forests of Asia and Africa (Sun et al., 2019). Hornbills are large birds belonging to the Bucerotidae family and out of the 59 extant species of hornbills, 31 are found in Asia (Franco and Minggu, 2019). Hornbills are often called as farmers of the forest in terms of their habitat, food, nest-site, and seed dispersers [7]. They are known as good indicators of the health of forests [8]. The unique breeding biology of these birds means that they are dependent on big trees of primary forest [9].

Bhutan provides a safe home for four hornbill species: wreathed hornbill Aceros undulates, oriental pied hornbill Anthracoceros albirostris, rufous-necked hornbill Aceros nipalensis, and great hornbill Buceros bicornis [10]. Asian hornbills are hunted for their body parts (casque and tail feathers for traditional attire), for the consumption of their meat, and for their body fat, which is believed to have medicinal properties [11, 12]. The population of hornbills is declining in an alarming rate due to continue habitat fragmentation and degradation, deforestation, logging, and increase in hunting pressure [13]. Today, their natural habitat remains largely fragmented and with an extraordinarily high level of threats to their persistence [14]. The target species rufous-necked hornbill has been reported to be extinct in Nepal and close to extinction in Vietnam [15, 16]. To date, no study has been carried out to know regarding hornbill’s habitat suitability and distribution in Bhutan. Against the backdrop, this study is intended to fill the existing gap with the objective to delineate the potential habitat suitability and identify critical variables associated with rufous-necked hornbill and great hornbill in Bhutan.

2. Material and Methods

2.1. Study Area

Bhutan is a small landlocked country spreading over an area of 38,394 Sq. km (Figure 1) and extending between latitudes 27° 31′ 53.11″ N and longitudes 90° 26′ 9.07″ E in the Eastern Himalayas which is bordered by China and India. The annual temperature ranges between 10°C and 24°C and annual precipitation between 300 m and 6000 m across Bhutan [17]. The territory of Bhutan is the home for 5500 species of vascular plants, 46 varieties of rhododendrons, 400 type of lichens, 430 varieties of orchids, and 200 types of forest mushrooms [18]. The human population is 735,553 of which 62.2% lives in rural areas, and their livelihoods depend on agriculture and livestock farming [19]. Additionally, more than 95% of Bhutan remains vegetated, of which approximately 70% constitutes natural forest cover [20]. The bird species of the country stands at 748 species of which 31 are globally threatened and 18 are part of the 37 endemic bird species which makes stronghold of bird diversity [21]. The protected areas system of Bhutan comprises five national parks, four wildlife sanctuaries, and one strict nature reserve [22].

3. Methods

The methodology outline is shown in Figure 2. A total of 56 occurrences points of two sympatric hornbills, i.e., great hornbill (n = 31) and rufous-necked hornbill (n = 25) were collected from the field survey. The field survey was conducted for four months (December 2018–March 2019). During the field survey, hornbills were sighted, and GPS coordinates (latitude and longitude) were recorded along with photo proofs. After the field survey, we used 21 environmental variables (BIO_1 to BIO_19, mean actual evapotranspiration (AET), and global aridity index (AI) (https://www.worldclim.org) were downloaded from the WorldClim dataset at a resolution of 1 km2. Topograhic variables, namely, elevation, aspect, and slope were derived from the digital elevation model of Shuttle Radar Topographic Mission (SRTM) having the resolution of 30m.

A test of autocorrelation was done using DIVA-GIS to remove correlated sites. After the test, 51 presence locations are retained and 5 autocorrrelated points were discarded (2 presence location for great hornbill and 3 presence location for rufous-necked hornbill. Likewise, correlation test of 24 variable layers was conducted using the ENM Tools version 1.3 [23] to test multicollinearity between predictor variables. The criteria adopted to predict multicollinearity between the variables was Pearson’s correlation coefficients (r). Pearson’s correlation coefficient (r) is a measure of the strength of the association between the two variables. After running the correlation test, variables with correlation coefficient value higher than 0.75 were discarded from the further analysis. Bio_1, Bio_2, Bio_4, Bio_5, Bio_6, Bio_7, Bio_8, Bio_9, Bio_11, Bio_12, and Bio_19 were removed and rest 13 variables were selected generating the model. To generate the habitat suitability map, all the spatial data layers were resampled at  m spatial resolution using the nearest neighbor re-sampling technique using ERDAS 2015 and ArcGIS 10.5 software. Then, the layers were converted to ASCII layer to feed into the MaxEnt 3.4.1. The model was iterated 100 times with 20 replication using the bootstrap to avoid any biases or extreme value that may arise and output was achieved as the mean among them. All rest of the settings we kept as default in the MaxEnt. The MaxEnt output of ASCII format are converted to raster and later categorized in ArcGIS for further evaluation.

The suitability map obtained will have the value between 0 and 1 where 0 represents least suitable and 1 represents most suitable habitat for the species. The maps were further classified into 4 categories: most suitable (0.6–1), moderately suitable (0.4–0.6), marginally suitable (0.2–0.4), and least suitable (0–0.2) (Yang et al., 2013).

4. Results

4.1. Potential Suitable Habitat Analysis of Great Hornbill (GH) and Rufous-Necked Hornbill (RNH)

The habitat suitability map of GH and RNH is divided into four classes as highly suitable, moderately suitable, marginally suitable, and least suitable. Each following classes are indicated by red, yellow, blue, and green, respectively. The MaxEnt result showed that the area wise class percentage analysis for GH comprised of 2%, 8%, 10%, and 80% (Figure 3; i, ii) whereas RNH comprised 3%, 6%, 15%, and 76% (Figure 3; iv, v) representing highly suitable, moderately suitable, marginally suitable and least suitable respectively. The predictor variables performance with the highest percentage contribution and permutation importance in predicting the species presence location for GH were elevation (30.4%), precipitation of wettest month (Bio_13) (27.2%), and aridity index (19%) (Figure 3; iii). In case of RNH, precipitation of wettest month (Bio_13) with 44.7%, aspect (22.3%), and aridity index (11.5%) (Figure 3; vi) were the highest contributor. The also model predicted that the highly suitable habitat for GH is located in the southern and central region and rufous-necked hornbill is high in northwest and north aspects of the country (Figures 4 and 5).

The curves below in Figure 4 illustrates how the logistic prediction changes as each environmental variable are varied while keeping all other environmental variables at their average sample value. Red lines are the mean response of the 100 replicates and blue is the ± one standard deviation. The x-axis represents the ranges of values of the environmental variables whereas y-axis represents the probability of occurrences on the scale (low probability–high probability). The graphs represent the effect of an individual environmental variable on the distribution of the great hornbill (RH) and rufous-necked hornbill (RNH). Precipitation of wettest month identifies the total precipitation that prevails during the wettest month in mm. The probability of occurrence of GH and RNH is high in 800–850 mm and 450–850 mm (Figures 6(a) and 6(b)). The range of BIO_13 in the study area is from 79–1181 mm. The probability of occurrence of both GH and RNH is high in slightly higher precipitation of the wettest month. The contribution of the aridity index in predicting the species presence location for GH and RNH was 19% and 11.5%, respectively (Figures 6(c) and 6(d)).

The response curve shows that the probability occurrence of GH is high at 14000 AI (Figure 6(c)), it shows GH prefers humid climatic condition. However, the response curve of RNH shows negative relationship with the aridity index (Figure 6(d)). The probability of occurrences of RNH is high in semiarid and arid region.

The contribution of elevation in predicting the species presence location of GH was highest (30.4%). The curve shows a negative relationship, i.e., elevation increases, the probability of occurrences of species decrease (Figure 6(e)), and we can see species occurrence probability is high at the elevation below 3000 m. In case of RNH, the contribution of aspect in predicting the habitat suitability was 22.3%. The aspect ranges from 1–360 degrees representing the direction between north, east, south, and west. The probability of occurrence of RNH is high in 340–360 degree (Figure 6(f)). It shows that the probability of occurrence of RNH is higher in north-facing aspect or sun direction. The jackknife test of variable importance shows the highest gain when the variable elevation is used in isolation (Figures 6(g) and 6(h)). The most influential environmental variable that decreases the gain is aspect. This indicates aspect is more effective in defining the model than any other variable.

Finally, the testing average AUC for 100 replicates runs of the model was AUC = 0.915 with the standard deviation 0.023 for the great hornbill and incase of rufous-necked hornbill was AUC = 0.872 with the standard deviation 0.060 (Figure 7).

5. Discussion

Birds are thought to perform the widest range of ecological functions in the ecosystem [24]. Birds benefit humans by providing a variety of important ecosystem services, including provisioning (food, feathers), regulating (seed dispersal, pollination), cultural (art, religion), and supporting (soil formation, nutrient cycling) services [25]. Hornbills are frugivore, which means they eat fruit. Frugivores play a critical role in the ecosystem’s structure and function [26]. However, due to the gradually increasing human population following anthropogenic disturbance, there has been a severe decline in the population of these birds.

Species distribution models (SDMs) are widely used in biogeography, biodiversity, and macroecology research to model species’ geographic distribution based on the correlation between known occurrence records and environmental conditions at occurrence locations [27]. SDM creates geographic maps of environmental suitability for a species [28]. We chose the MaxEnt model to predict the distribution of the great hornbill and the rufous-necked hornbill because its distribution prediction provides a powerful new tool that uses only presence data for species distribution modeling. Since its inception, MaxEnt has grown in popularity [29]. MaxEnt is a presence-only model that allows scientists to tap into the wealth of natural history data while avoiding the high cost of sampling species across their entire range [30]. Data on presence are plentiful, but data on absence are difficult to come by and frequently unreliable due to a lack of survey effort [31]. MaxEnt compares the distribution of presences along environmental gradients to the distribution of background points, which are drawn at random from the study area [29]. Furthermore, it considers interactions between environmental variables and appears to perform reasonably well with small amounts of occurrence data when valid occurrence data and appropriate predictor variables are chosen [30].

The distribution model for the great hornbill and rufous-necked hornbill was created using 21 environmental variables and three topographical variables (elevation, aspect, and slope). However, after performing a multicollinearity test in the ENM tool, only thirteen variables were used to prepare a habitat suitability model and conduct further analysis. The results revealed that the most important environmental variables for predicting the species presence location for GH and RNH were elevation (30.4%), precipitation of wettest month (27.2%), and aridity index (19%), and precipitation of wettest month (44.7%), aspect (22.3%), and aridity index (11.5%), respectively.

The probability of GH occurrence and elevation has a negative relationship, meaning that the likelihood of seeing great hornbills decreases as elevation rises. According to the current study, the probability of GH occurrence is high below 3000 masl. The great hornbill is known to frequent wet evergreen and mixed deciduous forests, ranging out into open deciduous areas to visit fruits and ascend slopes to at least 1,560 masl in south India [32] and up to 2,000 masl in Thailand [15]. The RNH is found in primary subtropical evergreen and deciduous forests between 600 and 2000 masl worldwide but have been seen as high as 2900 masl [33]. RNH presence was predicted in the species occurrences predicted map at Lhuentse, Yangtse, Punakha, and even up to Gasa, which is over 3000 masl. Habitat suitability analysis shows that GH and RNH are highly suitable for 2 and 3% of the total country’s area, respectively. Their occurrence and distribution are affected by factors such as vegetation type, habitat, temperature, availability of food, forest size, and tall old trees with holes [34].

6. Conclusion

The distribution potential of great hornbill and rufous-necked hornbill was delineated using species distribution modelling. This model is a very effective tool for mapping the suitable habitats of different species under the influence of topographic and climatic factors. In comparison, great hornbill has 2% of the area of Bhutan suitable for habitat whereas rufous-necked hornbill has total of 3% of most suitable habitat in Bhutan. However, rufous-necked hornbill habitat range is high in central, southeast, and northeast region of Bhutan. The probability of occurrence of the great hornbill is high in the southern region and elevation below 3000 msl. Rufous-necked hornbill prefers semiarid and arid region whereas great hornbill prefers humid condition. The model also revealed that the highest contributing variables for predicting this threatened species were aspect and precipitation of wettest month. The baseline information generated from the present study can be useful facilitating fieldwork, planning, and future scientific management of these species in the country.

6.1. Recommendation

Due to the limited time period, the presence data of the species was low and could not collect the distribution information in a different season. The future research on this species can be improved with data collection on the distribution of great hornbill and rufous-necked hornbill throughout its range of both breeding and nonbreeding season so that an appropriate conservation measure and better model could be undertaken.

Data Availability

The environmental data were obtained from the WorldClim dataset at a resolution of 1 km2 (https://www.worldclim.org). Topograhic variables, namely, elevation, aspect, and slope were derived from the digital elevation model of Shuttle Radar Topographic Mission (SRTM). The primary data that support the findings of this study are available from the corresponding author upon the request.


This work is made available as authorea preprint (Puri et al., 2022).

Conflicts of Interest

The authors declare no potential conflicts of interest.

Authors’ Contributions

NirKumar Puri, AmitKumar Verma, and Ramesh Chhetri performed the conceptualization, data curation, and software analysis. NirKumar Puri and AmitKumar Verma investigated the research work. NirKumar Puri wrote the manuscript. AmitKumar Verma, Ramesh Chhetri, Harish Bahadur Chand, and Sandip Rijal supervised the study. AmitKumar Verma performed the formal analysis. Ramesh Chhetri carried out the methodology and modified into paper. Ramesh Chhetri, Harish Bahadur Chand, and Sandip Rijal reviewed and edited the manuscipt.


The authors are very thankful to the Forest Research Institute, Dehradun, for creating such platform for them to learn, understand, and support for the completion of this work. The authors express their heartfelt gratitude to Ugyen Kelzang, Prakash Rai, and Wangchuk for their continuous support and help during the field survey.