Do Gibbons Live in Dense Forest? A Case Study from Upper Brahmaputra Valley, Assam, India

The present study aimed at predicting the potential habitat of Western Hoolock Gibbon (Hoolock hoolock) in the upper Brahmaputra river valley, Assam and to identify the conservation priority areas for systematic conservation of the species taking canopy cover into account. We used the maximum entropy algorithm for the prediction of the potential habitat of Gibbons using its current distribution with 19 environmental parameters as primary predictors. Spatio-temporal analyses of the habitat were carried out using satellite-based remote sensing and GIS techniques for two decades (1998–2018) along with Terra Modis Vegetation Continuous Field product to examine the land use land cover (LULC), habitat fragmentation, Normalized Difference Vegetation Index (NDVI) and tree cover percentage of the study area. To identify the conservation priority area, we applied a cost-effective decision-making analysis using systematic conservation prioritization in R programming. The model predicted an area of 6025.1 km 2 under high potential habitat, and a major part of it was found to overlap with the dense forest (80.04%) followed by moderately open forest (73.90%) and open forest (65.82%). The LULC change matrix showed deduction of forest areas in the predicted high potential habitat during the study period, while agricultural class showed an increasing trend. The fragmentation analysis indicated that the number of patches and patch density increased from 2008 to 2018 in the ‘very dense’ and ‘dense’ canopy regions of the gibbon habitat. Based on the conservation priority analysis, 640 km 2 area has been proposed to conserve a minimum 10% of the gibbons’ habitat. The current analysis revealed that in Upper Brahmaputra Valley most of the areas under dense forest and dense canopy have remained intact over the last two decades, at least within the high potential habitat of gibbons independent of the degree of area change in forest, agriculture and plantation.


Introduction
Primate diversity and abundance in a forest site are affected by the structural variables or quality of the habitat and indirect/direct anthropogenic impacts (Brown et al. 1985, Rylands 1987, Chapman and Peres 2001. Investigation of primate habitat and its insularisation is a major research concern in tropical forest ecosystem (Anderson et al. 2007), as almost 90 % of Primate taxa, in general, are facing highest threat of habitat fragmentation owing to its dependency on tropical forest (Andrén 1994, Marsh 2003. The species-speci c habitat requirement of a species is often considered to have correlation with habitat quality and feeding strategy of the concerned species. Therefore, it is important to understand the ecological exibility or limits of primate communities and species, in order to implement effective management strategies for their future conservation (Harcourt 1998, Lindenmayer 1999, Harcourt 2002). In the Indian context, unfortunately, primates are rarely studied in terms of habitat requirements.
Moreover, the major threats to the primates such as habitat loss, habitat fragmentation, and unsustainable economic development are unabated essentially obscuring the conservation processes.
Continuous monitoring of habitat and qualitative assessment of species' preferred habitat parameters is an urgent need to improve prioritizing conservation action. In this paper, we present a case report of Western Hoolock Gibbon (Hoolock hoolock) in the Upper Brahmaputra valley with special reference to habitat monitoring.
The species of the genus Hoolock are found in different forest types of northeast India (Das et al. 2003(Das et al. , 2009), Bangladesh Feeroz 1992, Das et al. 2003), Myanmar (Brockelman et al. 2009), and South China (Fan et al. 2011). Although the habitat characteristics of hoolock gibbons are strongly associated with the abundance of food, it is recognized that they prefer dense canopy forest areas in the tropical forest of northeast India in particular. Due to the high rate of forest fragmentation and habitat degradation in gibbon habitat, the loss of dense canopy is prominent posing a threat to the survival of Western Hoolock Gibbon (Hoolock hoolock) present in the south bank of the Brahmaputra River, Assam. Studies on gibbon in India are mainly concentrated on the species' population and behavioural aspects with very little emphasis on the quanti cation and monitoring of habitat parameters. The frugivory in gibbons are well established (Leighton 1987, Palombit 1992) with known food niche throughout its distribution range. However, minimal efforts have been realized in the associated forest structure evaluation and monitoring. Some research articles and PhD theses are available for gibbon population status and behavioural ecology particularly from northeast India that provide baseline ecological data of the species (e.g., Kakati 1997, Das 2002, Kakati 2004. As a number of primate species coexist in the region (Sharma et al. 2013), it is very important to determine the speci c habitat requirement of the species. It is noteworthy that studying sub-populations of gibbons in various pockets of northeast India's forest cannot be conclusive in many ways. Therefore, a landscape level study on habitat with the application of advanced tools and techniques is the need of the hour. Currently, no studies have been reported on gibbon habitat monitoring considering temporal satellite data from the region.
The present study aimed to investigate the site-speci c preference of gibbon habitats considering the canopy cover and forest cover by integrating maximum entropy algorithm to predict potential habitat and geospatial tools to incorporate major habitat forms. The study also emphasized to test the hypothesis that hoolock gibbons generally prefer dense forest or closed canopy. This information will provide quanti ed baseline data for the future conservation of gibbon in the Upper Brahmaputra River stretch in Assam, India. Subsequently, we identi ed fragmented patches in the gibbon habitat that need to be prioritized for the future conservation of the species. We adopted integer linear programming techniques (prioritizr R package) for the decision-making analysis, considering the canopy tree as a cost for the conservation (Hanson et al. 2018).

Study Area
The study was conducted in the Upper Brahmaputra valley as the hoolock gibbon is distributed on the south bank of the mighty Brahmaputra River. The southern part of the Upper Brahmaputra Valley comprises mainly four districts of Assam viz., Jorhat, Sivsagar, Dibrugarh, and Tinsukia. The study area lies between 94 0 30' E to 96 0 0'E longitudes and 27 0 0'N to 27 0 45'N latitudes covering an area of 9851.20 km 2 (Fig. 1). Four major protected areas including a National Park fall under this area which signi es the importance of the landscape in terms of oral and faunal conservation priorities. The vegetation of the study area is identi ed as wet evergreen forests of the Dipterocarpus-Mesua series (Champion and Seth 1968). The climate of the Brahmaputra valley as a whole is similar to the Southeast Asiatic Monsoon climate modi ed as per local physical conditions (Deka et al. 2013). The valley receives as high as 3900 mm rainfall annually in the extreme northwest and northeast hilly tracts (Nath and Deka 2010). The mean annual temperature ranges from 23 0 C to 24 0 C.

Methods
The survey was conducted from January, 2017 to December, 2018 for gibbon occurrence and ground veri cation (ground-truthing). Gibbon locations of both direct sightings and heard calls were collected through handheld GPS. Along with the occurrence locations, the land use pattern was also noted down to prepare LULC supervised classi cation signature. Occurrence data collected were used in the predictive model thereafter.

Predictive Distribution Model
Maxent software version 3.1.0 (Computer Sciences Department -Princeton University, 2004) was used for predicting the potential distribution area of H. hoolock. Nineteen Environmental layers comprising temperature, precipitation, etc. (Supporting information; Appendix S1) and 23 species occurrence data were fed into Maxent, and the nal predicted potential distribution area was projected using ArcGIS 10.3.
The detailed explanation of ENM methodology using Maxent had been described by Phillips et al. (2006).
The occurrence data of H. hoolock was gathered through extensive eld surveys conducted along the Upper Brahmaputra valley, Assam, India. A total of 23 independent distribution localities of gibbons were collected from the eld, and all localities were used in the nal modelling process. Nineteen environmental layers including temperature, precipitation, etc. were downloaded from the WorldClim website (www.worldclim.org) with spatial resolution of 30´´(seconds).
The Jackknife validation methodology was performed following the method developed by Pearson et al. (2007), which was shown to be effective for sample sizes of 25 or less. To avoid over tting of the test data, the regularization multiplier value was set at 0.1 (Phillips et al. 2004). The maximum number of background points were 1000. Linear quadratic and hinge features were used. We selected 80% of the data for training and the rest 20% for testing and 10 percentile threshold rules was employed. A total of 100 runs were set for model building (Flory et al. 2012). AUC (Area Under the receiving operator Curve) was used to test the model's goodness-of-t. The highest AUC was considered as the best performer. The contributions of the variables were assessed through the Jackknife procedure. The nal output was divided into three potential distribution areas that were regrouped with a range of 0-1: Low potential categories were decided viz., agriculture, grassland, forest, plantation, sandbar, settlement, and water (Supporting information; Appendix S3). Vector polygon was drawn around each representative predetermined LULC categories (collected from ground-truthing) to prepare the spectral signatures with minimal confusion points (Gao and Liu 2010). A total of 70 spectral signatures for the respective LULC categories derived from the satellite imagery were recorded using the pixels enclosed by these polygons.
These spectral signatures were then used to reclassify the images using a maximum likelihood classi cation algorithm that classi es the pixels based on the maximum probability of belonging to a particular class (Richards and Jia 2006). Then, the classi ed images were ltered using a neighbourhood majority function.
The LULC change-transition matrix was computed using the overlay procedure in ArcGIS to quantify the area converted from a particular LULC class to another LULC category during the study period. Later, visual interpretation was used to address the mixed pixels problem (Twisa and Buchroithner 2019).

Accuracy assessment
The non-parametric Kappa test was performed to measure the extent of the classi cation accuracy to account for diagonal elements and elements in the confusion matrix (Rosen eld and Fitzpatrick-Lins 1986). A confusion matrix was constructed with each row representing LULC classes in the classi ed map and each column representing the reference LULC classes. The kappa co-e cient was determined from this matrix for each classi ed map. Kappa Co-e cient is the degree of agreement or precision between the classi ed map and the reference data (Rosen eld and Fitzpatrick-Lins 1986, Congalton 1991).

Percentage Canopy Cover
Terra Modis Vegetation Continuous Field (VCF) product was used to extract the percentage tree cover of the study area. The VCF products is a monthly composite of Terra MODIS 250m and 500m Land Surface Re ectance data, and the product provide a gradation of tree percent cover (percent of a pixel covered by canopy), percent non-tree vegetation (non-tree canopy pixel), and percent non-vegetated (pixel with no vegetation) (https://lpdaac.usgs.gov). We constructed a composite image of tree canopy cover (mean aggregated) for the year 2008 and 2018 on the Google Earth Engine platform (GEE). GEE is a platform that allows users to compute its multi-petabyte catalogue of satellite imagery and geospatial dataset with planetary-scale, and export analysis in Geotiff format or tabulated data (Gorelick et al. 2017). We exported the percentage tree cover (250 m spatial resolution) as a TIFF image and reclassi ed it into three categories following FSI forest cover classi cation (www.fsi.nic.in) viz., no canopy (0-10%), open canopy (10-40%), dense canopy (40-70%), and Very dense canopy (>70%). We used ArcGIS for the reclassi cation and extraction of pixel data. To minimize the over prediction of Maxent model in terms of conservation planning, prioritizR (Hanson et al. 2020) package in R platform was used to nd out the potential feasible areas for conservation action within the high potential zone of the predicted distribution map. A 4 km x 4 km grid map was overlaid on the high potential distribution zone and grids that were in protected areas were excluded using the 'locked out' function in the conservation problem de nition. Also, grids with more than 30% human settlement areas were excluded for feasible conservation planning of gibbon protection in the region. The cost of conservation was assessed based on percentage canopy present in each grid in a scale of 1 to 5 (5=0-20%; 4= 20-40%; 3=40-60%; 2= 60-80%; 1=>80%). The lowest canopy percentage grid was assigned a scale of 5 and the highest was assigned 1. Three different threshold percentages accounting 10%, 15%, and 20% species distribution area were intended to conserve and the problem was de ned accordingly.

Patch analysis and conservation prioritization
The output maps were then checked for irreplaceability and nal maps were prepared. The output predictive map has shown that the high potential habitat zone comprised the highest area (6025.1 km 2 ) followed by moderate (2003.47 km 2 ) and low potential (1729.19 km 2 ) habitat categories ( Fig. 2 [b]).

LULC Change Matrix
Accuracy assessment showed that the kappa values of the LULC classi cations were 0.764, 0.693, and 0.757 for 1998, 2008 and 2018 respectively. LULC analyses revealed that agriculture, grassland, and settlement classes covered majority of the study area (9851.2 km 2 ; Table 1; Appendix S4). Plantation and settlement classes showed an increase in area percentage between 1998 and 2018 indicating a signi cant role of human in land conversion. Urbanization and population expansion are major reasons for the continual surge in settlement areas (102%; Table1). On the contrary, forest areas showed a slight decline (4.5%) during the two decades. Area decrease in forests is expected to be higher, but due to the inclusion of scrublands in the forest class in the current study, the resultant area decline might have shown lower.  The distribution of forest cover among these thresholds was analysed and it was observed that the highest percentage of forest cover was under high potential habitat zone in the year 1998 which gradually decreased in 2008 and 2018 (Fig.3a). Another important land use category i.e. agriculture was highest in 2018 under high potential habitat area (22.71%) of gibbon, that has shown an increasing trend from 1998 onwards (Fig.3b). A similar trend was followed by moderate potential (6.41%, 20.89%. 21.41% for the years 1998, 2008, and 2018 respectively). However, in the low potential habitat, agriculture covered the highest area percentage in the year 2008 (26.58%). The plantation area displayed a minimal decrease in area from 1998 to 2018 in all the three potential thresholds (Fig. 3c).

NDVI and Potential Gibbon Distribution
The nal NDVI map of 2018 (Supporting information; Appendix S7) has shown that, the maximum tea plantation area (1148.70 km 2 ) has been included in the total dense forest area (4312.14 km 2 ) that is classi ed from NDVI. Lesser areas of tea plantation were included in moderate (58.92%) and open forest (60.65%) (Supporting information; Appendix S8).
Overall areas of NDVI classes under various potential habitat areas of gibbon in the study area revealed that the highest percentage area (80.04%) of dense forest was under high potential habitat (Fig.3d). In moderate potential threshold, 73.90% dense forest was included, and 65.82% dense forest was found in low potential threshold. While comparing NDVI classes omitting the areas of tea plantation in predictive distribution thresholds, it was found that 76.21% dense forest area under high potential threshold (Fig.  3e). In moderate potential threshold, 70.83% dense forest was recorded and 63.69% dense forest was found in low potential threshold. Thus, the dense forest with tea plantation was higher in percentage in high potential habitat zone (Supporting information; Appendix S9).

Protected Area Coverage
Overall, the protected area covers 20.82% of the study area (Supporting information; Appendix S11). The protected area covers 17.73% of the total high potential gibbon habitat whereas moderate and low potential habitats of gibbon in the study area include 27.51% and 23.74% of protected forest respectively.
Out of the total, very dense canopy in the study area 95.03% area was within protected range in 2018 whereas the same was 97.19 % for the year 2008. But when compared to the percentage of the total area irrespective of canopy categories only 7.62% very dense canopy was found in 2018 within the protected areas. Likewise, for the year 2008, the percentage of the same was 6.27% (Supporting information; Appendix S12). Furthermore, within the protected range, 75.35 km 2 has been lost from dense canopy from 2008 to 2018 whereas the very dense canopy increased by 133.17 km 2 during the same decade. On the contrary, in the non-protected range, a signi cant increase in dense canopy forest was recorded (548.25 km 2 ) from 2008 to 2018. The very dense canopy was also found to increase in the non-protected range (21.43 km 2 ), but the area addition was comparatively lower than that of the corresponding protected range.

Patch analysis and conservation recommendations
The FRAGSTATS result showed an increase in the number of fragments in both dense and very dense canopy forests from 2008 to 2018 within the high potential gibbon distribution area (Supporting information; Appendix S13). In 2008, the mean patch sizes in dense and very dense canopy areas were The proposed conservation planning to conserve 10% of the high potential habitat excluding both existing protected areas and high human habitation areas (>30%) showed 40 grids (640 km 2 ) highlighted for conservation planning (Fig. 4 A-F). Similarly, 58 grids (928 km 2 ) and 81 grids (1296 km 2 ) were highlighted to protect 15% to 20% of species' habitat respectively.

Discussion
Gibbon habitat and monitoring is crucial in tropical forests owing to its declining population trend throughout distribution range. The present study revealed the potential impact of forest dynamics in gibbon distribution and survival in the Upper Brahmaputra valley located along the south bank of the Brahmaputra River in Assam. The study area is one of the major habitat zones of western hoolock gibbon in India. Although most of the forested areas of this zone are under legal protection, the heterogeneity of the forest composition and anthropogenic intrusion render the agreement notably. The LULC analyses revealed that Agriculture, Grassland, and Settlement classes covered the majority of the study area. Plantation and Settlement classes showed an increase in area percentage between 1998 and 2018 suggesting a substantial anthropogenic role in land conversion. Urbanization and population expansion are major reasons for the continual surge in settlement areas. On the contrary, forest areas showed a slight decline (4.5%) during the two decades. Area decrease in forests is expected to be higher, but due to the inclusion of scrublands in the forest class in the current study, the resultant area decline might have shown lower. It is also reported that fragmentation is very severe in the lowland rainforest of the Upper Brahmaputra valley and one-third of the forest cover area has been lost in the last century (Sharma et al. 2013). The higher percentage change in forest cover between 1998 and 2008 depicted the ineffective efforts of habitat conservation of the species in the study area. As gibbons are known to prefer dense or pristine forest for living, it is a major concern for conservationist to protect dense forest areas. The shrinkage of forest cover areas in fragments of known gibbon habitats in the study area has also been reported in earlier studies (Kakati 2004, Sharma et al. 2012a, Sharma et al. 2013) due to which the population is eventually exposed to open forest. Therefore, in the present study, the predictive potential habitat of the species is used as a primary tool to compare the habitat of the species on a temporal scale with particular reference to forest covers and canopy covers. Looking at the landscape and other sociopolitical disturbances in the study area, remote sensing methods used to track biodiversity have also proved to be effective (Turner et al. 2003, Corbane et al. 2015. Nevertheless, ground veri cation is a necessity for better con dence in the results. The population data on western hoolock gibbon is readily available for the study area as the Upper Brahmaputra valley is considered as the hotspot of gibbon distribution in India (e.g., Kakati 1997, Das 2002, Kataki 2004, Kakati et al. 2009). However, the study area is signi cantly lacking information on major habitat parameters of the species.
The Maxent model predicted that the area under the high potential habitat zone in the study area is as high as expected. It is noted here that Maxent does not account for the historical distribution of species nor does it recognize physical barriers of the restricted-range species (Peterson et al. 1999, Soberón andPeterson 2005). This may be the explanation for the inclusion of new areas of the species' high potential habitat dependent on bioclimatic suitability. The same argument was also supported by Sarma et al. (2015) in . On the contrary, the present study also con rmed that the highest percentage area of Dense Forest was found under the high potential habitat zone of gibbons. It is noted here that tea plantation is the key driver that affects the NDVI map of the study area as tea plantations are classi ed under Dense Forest categories. Since we have omitted the tea plantation areas from Dense Forest zone by testing the output map with ground data, the high potential habitat zone still represents the maximum Dense Forest area even though the tea plantation area is removed.  (Chetry et al. 2007), the present study area includes as many as 42 RFs. This study therefore strongly supports the implementation of stringent laws on small fragments of rainforests in the Upper Brahmaputra Valley for the potential survival of gibbons, taking into consideration that the species prefers dense forests. It should be mentioned here that gibbons are still present in small isolated fragments of the study area, irrespective of the canopy structure or the condition of the habitat. In view of the narrow geographical range of western hoolock gibbon, the species must survive with a reduced abundance of resources before ecosystem conditions are adequately controlled.
Do gibbons live in dense forest? Our analysis shows that over the last two decades in Upper Brahmaputra Valley, Assam, most of the areas of dense forest, as well as dense canopy, have remained intact, at least within the high potential habitat of gibbons independent of the degree of change in forest, agriculture and plantation covers. The protected area network within the study area is of greater value in ensuring improved protection of the remaining gibbon habitat in the study area. However, the investigation of other possible causes, such as food availability and food preference, is often a critical concept for recognizing the species' preferred habitat. The decadal changes in forest cover and canopy structure of the study area therefore may be a few of the many drivers that affect habitat preferences of gibbon. On the other hand, knowledge of the population dynamics of the species in the study area can also lead to recognizing the need to study unique viewpoints of the species' ecology as the species is included in the endangered category of the IUCN Red List. Here, in this study, we did not consider the physical conditions of the habitat and the constant land-use pressure posed by the human population on the habitat of the species.
Though proposing future conservation areas within the predicted high potential zone of gibbons, we deliberately omitted the current protected areas believing that the species within protected areas are still protected under forest preservation legislation. The human settlement areas have already been ltered to ensure ease of implementing the Conservation Action Plan. However, the cost of the conservation was objectively identi ed based on the fact that gibbons prefer dense canopy for survival. We also believe that the outcome of this study would be the basis for the comprehensive investigation of factors in uencing the responses of gibbons to environmental and climatic changes across their distribution range. Declarations