Changing landscape configuration demands ecological planning: Retrospect and prospect for megaherbivores of North Bengal

Summary

The Gorumara National Park (GNP) is an important conservation area located in the northern region of West Bengal State, India, as it provides habitat for three megaherbivores: Indian One-horned rhinoceros (Rhinoceros unicornis), Asian elephants (Elephas maximus) and Gaurs (Bos gaurus). It harbours one of the last population of the one-horned rhino. In the present study, landscape change and configuration were investigated by comparing three Landsat images, from 1998, 2008 and 2018. The images were classified into six different landcover classes following standard methodology. The present study also involves evaluation of landscape and anthropogenic predictors influence on the megaherbivores of GNP, followed by future landcover simulation for the year 2028. The result shows a significant decrease in the grassland cover from 18.87 km² to 8.27 km² from 1998 to 2018, whereas the woodland cover has increased from 50.14 km² to 62.09 km² between 1998 and 2018. The landscape configuration indices such as Number of Patches (NP), Patch Density (PD), Interspersion and Juxtaposition (IJI), Aggregation Index (AI) and Mean Shape Index (SHAPE AM) indicated that the landscapes has lost complexity in the spatial placement of patches of different Land Use and Land Cover (LULC) classes. Also, the landscape over the three decades has become uniform in terms of diversity of patches, because of earlier plantation activities by the forest managers. Result also indicated that grassland, along with its class metrics are the top predictors contributing 43.6% in explaining the spatial distribution of megaherbivores in GNP. Results from the simulated landcover of 2028 suggest a possible decline in overall grassland by 6.23% and a subsequent upsurge in woodland by 6.09% from 2018. The present result will be useful in guiding the forest management in developing habitat improvement strategies for the long- term viability of megaherbivore populations of rhino, gaur and elephant in the GNP.

Keywords:

Grasses – Forests – Habitats – Conservation science – Wildlife – Grasslands – Surface water – Elephants

Introduction

Detecting the rate of change and pattern of change in the landscape has been considered as an essential theme in ecological research [1, 2]. Multi-temporal change analysis of remotely sensed data enables us to understand the dynamicity in landscape composition [3,4,5,6]. Understanding the pattern of change in the landscape configuration of Protected Areas (PAs) including National Parks (NP), Wildlife Sanctuaries (WLS) and reserves would eventually help the managers making strategies for effective conservation and management. A number of studies are available which have demonstrated the applicability of mapping and landscape change analysis for long-term conservation and management [7,8,9,10,11].

Megaherbivores (mammals with bodyweight over 10³ kg) generally need more food and space than mesoherbivores, and these two, being the common constraints in areas with high population pressure as there population regulation through predation processes is very limited [12]. In today’s world, most productive habitats like floodplains and foothills are the areas where the megaherbivores are generally found. However, such regions are profoundly affected by habitat loss, fragmentation and other anthropogenic disturbances, restricting their population into small pockets [12,13]. Moreover, megaherbivores are known to have powerful effects on both the assembly and composition of vegetation and hence play an essential role in ecosystem dynamics [14].

Gorumara National Park (GNP) is known for its large assemblage of megaherbivores i.e. Indian One-horned rhinoceros (Rhinoceros unicornis), Asian elephants (Elephas maximus) and Gaurs (Bos gaurus). The habitat stretch of (GNP) includes territorial forests of Jalpaiguri, Kalimpong and Cooch Behar Forest Divisions also Jaldapara Wildlife Sanctuary, Chapramari Wildlife Sanctuary. But the total stretch is not continuous, recently due to enormous expansion of tourism, roads and railways with other anthropogenic activities resulted in substantial degradation and fragmentation of the habitat stretch [15]. These activities may lead to a loss in both movement corridors as well as genetic vigour of the wildlife populations [16,17].

Despite various threats in the past decades, the park management could be able to reduce the poaching threat successfully, that has resulted in a considerable increase in rhino population, i.e. 50 individuals [18]. However, due to massive population expansion of other sympatric megaherbivores, mostly Bos gaurus with over 900 individuals and Elephas maximus in GNP have resulted in an ecological crisis and have led to increased interspecific competitions [15]. For this reason, rhinos are frequently found outside its natural habitat (far-flung as 50 km) to meet their foraging needs, and some of them are staying seasonally or even permanently [19].

Importance of grasses and grassland, both on their qualitative and quantitative aspects is deeply associated with megaherbivores of GNP. For Rhinos, they were found to prefer grasslands intermixed with wetland and riverine forests [20]. Studies also indicate that grasses are the most preferred diet for rhino [21,22,23]. Similar types of findings were also observed across most Indian one-horned rhino reserves from India and Nepal [24,25]. Hence it can be opined that the rhinos are more habitat specialised than other megaherbivores and depend on food quality rather than quantity [22,26].Moreover, in Indian subcontinent gaurs are primarily grazers, and their diet is largely composed of grasses. The preference of grasses by gaur is supported by findings from various studies throughout Indian subcontinent [27–31]. On the other hand, elephants are known to exploit a large variety of species, yet the dominance of grasses in the diet of the elephant was documented by various studies, including GNP [32–35].

However, among the many ecosystems of GNP, grassland are the most affected ones in recent times, which is evident by a sharp decline in both productivity and species diversity. Leading causes of which include overgrazing, soil erosion, nutrient depletion, salinization, pollution, disruption of hydrological systems, and conversion of natural areas into croplands, monoculture plantations and ill-planned developmental activities resulting in exceeded carrying capacity for rhino [19]. The carrying capacity for rhino in Gorumara was first assessed to be 20 in the year 1995 [36]. Followed by several other studies which suggest they exceeded carrying capacity threshold for rhino in GNP [19,37].

Landscape-level characterization remains associated with a qualitative and quantitative evaluation which helps in identifying the linkage between spatial pattern and ecological processes, like composition, configuration and connectivity [16,38,39]. Landscape composition mainly refers to the number of different elements, their spatial distribution, relative abundance and diversity, while landscape configuration refers to physical orientation, heterogeneity, and distribution within the landscape [39]. Globally, forested habitats, including forests in protected areas, are facing severe management challenges posed by various planned and unplanned degradation drivers [40,41].

Spatial heterogeneity of these landscapes influences the distribution, abundance and population dynamics of species [42]. As a consequence, new ecological theories along with new applications for planning and monitoring have been developed. Studies have also been conducted towards identification and evaluation of drivers which greatly influencing the spatial distribution and habitat suitability of the species [43,44,45]. Moreover, several researchers have also used landscape metrics and presented their imperativeness in developing a better understanding of the landscapes which species are utilizing [45,46,47]. Most prominently indices about patch structure, their spatial placement and configuration in the landscape are useful in understanding the ecology of large mammalian species including elephants, rhinos [48–51].

Furthermore, recent techniques in remote sensing have provided opportunities for developing models for further understanding of the future landcover dynamics, which will eventually help in making management decisions [52]. Among various approaches for landcover modelling and predicting land use in future scenarios, Markov chain analysis has been used most commonly [53]. Markov chain analysis is based on a stochastic approach which predicts a specific state of an extensive system at a particular temporal scale based on the state of past trends by using different topographic and landcover drivers [54,55].

Hence, the present study has been designed by combining Land Use and Land Cover (LULC) analysis of georeferenced satellite data of three decades (1998, 2008 and 2018) with landscape ecology theories and methods. Our overall goal was to generate a better understanding of LULC change patterns and configuration of the GNP landscape, which is facing tremendous pressure because of recent human colonization and related development. We have also evaluated the influence of landscape and anthropogenic predictors on the spatial distribution of megaherbivores in GNP. Further, to understand the future scenario, we have used Markov-cellular automata (CA) technique based Artificial Neural Network model (ANN) for Spatio-temporal modelling of future landcover for the year 2028.

We seek to develop baselines for a variety of landscape indices, which can describe the conversion of formerly contiguous habitats and its impacts on animals with relatively large territories including rhinos, elephants and gaur [56–58]. The Indian one-horned rhino is the flagship and also the top conservation priority species of GNP, and most of the management strategies of the park are mostly concentrated around this species. Our analysis of three decades of landscape change will provide a retrospective view on landscape change and will help in planning long-term conservation measures for the megaherbivore populations.

Material and methods

Study area

The GNP is located in the northern region of the West Bengal state of India (Fig 1), falling in the 7B-Lower Gangetic plan biogeographic zone situated between latitude 26°47’12.5” North to 26°43’25.6” North and longitude 88°52’04.2” East to 88°47’07.3” (Fig 1), with an area of 90 km² [59]. The park is located in the foothills of the eastern Himalayas, popularly known as Dooars, composed of flood plains, hills, forests, grasslands and plantations. The landscape is known for its rich biodiversity which includes 193 species of birds, 22 species of reptiles, 40 species of fishes and other macro and micro fauna [60,61]. The rich diversity of herbivores is represented by Indian Rhinoceros, Gaur, Asian Elephant, Chital, Sambar, Barking deer, Hog deer and Wild boar. Common Leopard is the only large carnivore residing in the study area.

The terrain of the landscape is differentiated into distinct plateaus and lacks close contours. The river system includes three main seasonal rivers, namely, the Murti, Indong and Garati Rivers, which drain into the Jaldakha River. The vegetation in the landscape is dominated by the tree species Shorea robusta, Lagerstroemia speciosa, Senegalia catechu and Dalbergia sissoo, and classified as riverine forests, sal forests, wet mixed forests and sal savannah [62]. The GNP contains a mosaic of mixed land-use and land-cover types in which human-wildlife conflict management is a major challenge for the authorities as well a source of threat to the large mammals in the study area [63].

Earlier forested habitat of GNP was a part of Jalpaiguri Forest division comprising of an area only 8.62 km² and the area was declared as Gorumara Game Sanctuary vide Govt. Notification No. 5181-For, Dt: 02.08.1949 just after independence. Afterwards through subsequent notification finally it got the status of National Park in the year 1998 [62,63]. Before 1994, the primary aim of this forest areas was to maximise timber production through the plantation of commercial species, such as Sal (Shorea robusta), Teak (Tectona grandis) and Jarul (Lagerstroemia speciosa). Hence, the forest continues to be exploited under various working plans from 1892 onwards until the end of the 4^th working plan (1926–27 to 1945–46). However, after the notification of the area as National Park, the objective has changed from production forestry to protection, for the conservation of wildlife and biodiversity. The first Management Plan of GNP was prepared and approved by Ministry of Environment and Forests, Govt. of India in 1997–98 with the aim of habitat improvement with focus on flagship species rhinoceros and other wild animals. As a habitat improvement strategy, canopy opening was done till year 2000, and then it has been stopped due to ban imposed by the Supreme Court of India. Later in the year 2006, the ban was lifted by the apex court for limited activities, but no further canopy opening initiative has been documented in the recent management plan [62,63]. Presently, the management objectives in GNP include conservation of biodiversity with special thrust on the rhino, development of GNP as elephant reserve, community participation in forest management and promote ecologically sustainable tourism for nature education and awareness creation.

Study design and data collection

For understanding the LULC change in the study area, we used three sets of Landsat satellite data (1998, 2008 and 2018), downloaded from the online USGS postal (https://earthexplorer.usgs.gov). For the landcover-change analysis, we used two types of Landsat sensor data. For 1998 and 2008, we used Landsat- 5 TM (Thematic Mapper), and for 2018 we used Landsat- 8 OLI (Operational Land Imager). All images were scaled for a spatial resolution of 30 m. To address the change in sensors from Landsat-5 TM and Landsat-8 OLI sensors, we used DEM data and geo-referenced the respective Landsat images following the robust method suggested by [28,29]. All images were projected in WGS_1984_UTM_Zone_45N (Table 1). Atmospheric, along with DOS1, corrections have been done for all three images using the semi-automatic classification tool for QGIS [64] followed by ML (Maximum Likelihood) technique for classification. All three images were classified into six major LULC classes: Woodland, Grassland, Shrubland, Bareland, Riverbank and Water. For the present study we have applied similar classification approach for the classification of 1998 and 2008 images by overlying training data polygon of 2018 on both 1998 and 2008 images followed by selection and elimination process for training samples which showed change in the cover type [65–67]. A total of 137 polygons comprising 973 pixels were taken during the field survey and were used on all three decadal images (1998, 2008 and 2018). Finally, we assessed the accuracy of the classified image by estimating the error matrix along with the kappa statistics. Sample points for accuracy estimation were also collected during the field visits; at least 60 samples location were taken for each land cover type for the present land cover.

Furthermore, Google Earth images and Forest Management Inventory data have been used for the past decadal images [68,69]. We calculated transition probabilities to determine the probability of transition from one land-cover type to another [68]. All image related analysis, as well as accuracy assessments, were carried out using the semi-automatic Classification plugin of QGIS, ENVI 5.1 and ESRI ArcGIS 10. 6.

We used the FRAGSTATS 4.2 program [70] to calculate a variety of landscape indices (S1 Table) for characterizing changes in land cover and landscape configuration across the three-decade images [71,72]. All the maps including the LULC, landscape configuration indices such as diversity index, patch density and patch density richness were developed using the ArcGIS 10.6 software program. The LULC and the landscape configuration change were assessed by comparing the three decadal data and by estimating the transition probability matrices between 1998 and 2008; 2008 to 2018 of the respective multi-temporal images using a Markovian approach.

Evaluation of landscape and class predictors influence on the Megaherbivores of GNP

Predictor inspection was performed by using Generalized Additive Model (GAM) [73] in SAHM package of VistTrails environment [73,74]. For understanding the relative response to landscape configuration, topographic as well as anthropogenic predictors by the mega-herbivores of GNP the GAM was chosen for its adaptability in dealing with a non-linear relationship between response and explanatory variables [75]. The land cover classes and linear features (Road and Rail) variables were prepared by converting raster into vector followed by generating euclidian distance, whereas, landscape and class level metrics were prepared by using the moving window function of FRAGSTAT Ver. 4.2 (S4 Fig).

A total of n = 97 spatially independent GPS locations of direct sightings as well as indirect signs for all three mega-herbivores (Elephant, Rhino and Gaur) were collected by surveying line transects marked systematically in the landscape after stratification of habitat types in GNP between the year 2016 to 2018 (Fig 1). For evaluating relative responses of different variables, pseudoabsence points have been assessed by generating random points in ArcGIS 10.6 in the entire landscape. For all presence records the GPS location, habitat type, distance to water and distance to rail and road have also been recorded.

Landcover simulation for the year 2028

We used Artificial Neural Network (ANN) based model for generating a future land cover map for the year 2028 of GNP in MOLUSCE plugin of QGIS using the classified images of 2008 and 2018 [68]. Furthermore, for simulating the landscape cellular-automata method based on Monte Carlo algorithm was used [68,76,77,78]. A total of five uncorrelated variables with Pearson’s correlation threshold less than 0.8 were selected as predictor variables for the simulation, and all the rasters have been resampled in 30 m cell size. The topographic variables, i.e. Compound Topographic Index [79,80], linear aspect and roughness, were used as surface texture/configuration drivers. Moreover, temperature and moisture drivers, i.e. Heat Load Index [81] and Integrated Moisture Index [82] were kept as a proxy for climatic drivers. All raster images were prepared using Geomorphometric and Gradient Metrics Toolbox Ver. 2.0 [83]. The validation of the simulated landcover of 2028 was carried out using percentage correctness along with kappa statistics.

Results

The three Landsat images were classified into six dominant LULC (Woodland, Grassland, Shrubland, Bare land, Riverbank and Water) (Fig 2). The overall accuracy for 2018 was evaluated to be 88.61% with a kappa value of 0.863 (Table A in S2 Table). The overall accuracy of the image 1998 and 2008 was 89.97% and 81.66% and kappa value was found to be 0.88 and 0.78 respectively (Table B and C in S2 Table). The comparative analysis of the three images revealed that woodland was the dominant land cover and increased in all three successive years (1998, 2008, and 2018), (Table 2), in contrast, grassland areas declined over the time. The bareland cover has increased from 6.72 km² to about 9.86 km² between 1998 and 2018, primarily due to loss of water bodies and narrowing of river lines.

The transition analysis from 1998 to 2008 indicates that the grassland cover has been converted to woodland at the most (0.517 transition value), followed by bareland to scrubland (0.375 transition value) (Table 3) (S1 Fig) (S3 Table). A similar trend was also observed from 2008 to 2018, where grassland has been converted to woodland cover (0.496 transition value). Furthermore, during 2008–2018 decade river bank got converted to bare land (0.451 transition value) and water cover has also been converted to river bank cover (0.323 transition value) (Table 4) (S1 Fig) (S3 Table).

Landscape configuration change analysis

The change in landscape configuration indices from 1998 to 2008 and from 2008 to 2018 indicates a decrease in the number of patches (NP) and patch density (PD) (Table 5). A smaller change was observed in NP and PD from 2008 to 2018. The mean patch area (AREA MN) increased by 24% and by 2.04% from 1998 to 2008 to 2018 respectively. Both Shannon’s diversity (SHDI) and evenness (SHEI) of patches decreased. The Aggregation Index (AI) and Interspersion and Juxtaposition index (IJI) of patches increased by 5.64% and 14.67% respectively between 1998 to 2018. The Edge Density (ED) found to be following a decreasing trend starting from a reduction by 39.50% from 1998 to 2008 and by 16.26% from 2008 to 2018 suggesting a decrease in the heterogeneity of the overall landscape in multitemporal scale. Results show that about 886.14 ha of grassland have been converted into the woodland from 1998 to 2008, followed by conversion of 555.57 ha between 2008 and 2018 (Fig 3). Other major insertions have occurred between grassland to shrubland by 217.62 ha. and by 328.95 ha. followed by shrubland to bareland by 291.78 ha and by 266.4 ha. between 1998–2008 and 2008–2018, respectively (S3 Table) (S1 Fig).

LULC class-level configuration change analysis

The proportional area of the woodland patches in the landscape (PLAND) increased by 11.94% from 1998 to 2018 (Table 6). The number of patches (NP) along with patch density (PD) for woodland has decreased from 86 to 26 and 0.94 to 0.28 from 2008 to 2018, respectively. The Largest Patch Index (LPI- Woodland) increased by 4.91% from the years 1998 to 2018 (Table 6). Thus, the park has more woodland in a few larger patches. Area weighted Mean Shape Index (SHAPE AM-Woodland), was highest in 1998, and then fell by 81% in 2008 and to 78.23% by 2018, indicating a reduction in shape complexity of woodland patches (Table 6).

In the case of grassland cover, the PLAND value was highest in 1998, at 18.87%, followed by 12.31% in 2008 and 8.27% in 2018 (Table 6). NP declined from 717 patches in 1998 to 389 patches in 2018. The mean patch area (AREA MN-Grassland) declined by 36.64% in 2008 and then increased by 29.69% by 2018. The shrubland, cover has decreased from 1998 to 2018, indicated by decreases in PLAND value between 1998 and 2008 by 23.06%. However, followed by a consecutive increase again by 18.73% in the year 2018. The number of patches (NP Shrubland) has increased from 259 patches to 339 patches in 1998 to 2008 initially and then decrease significantly to 277 patches in 2018. The distribution of small and isolated patches from the complex cluster of the larger patches was measured by Euclidean nearest neighbour distance metric (ENN_AM). The most contrasting change has been observed in shrubland cover, it had increased from 75.39 in 1998 to 111.82 in 2018.

The response of megaherbivores to landscape and class variables

The present results indicate that Grassland, along with its class metrics are the top contributors as it explains about 43.6% deviation in the spatial distribution pattern of megaherbivores of GNP (Fig 4). Out of which Area weighted mean (AREA_AM_Grassland) for grassland contributed highest followed by the euclidian distances to grassland. Moreover, grassland, along with its metrics shows a positive response in explaining the spatial distribution of megaherbivores (S3 Fig). Further, the landscape configuration variables of grassland cover including PD, NP, LPI and Area AM responded positively to the distribution of megaherbivores in GNP (S2 Fig) (S3 Fig). The Area AM of Shrubland cover contributes about 5.6%, suggesting the importance of shrubland for the megaherbivores in GNP. Among the anthropogenic variables, rail (railway line) contributes most by 4.9% and shows a negative response with the distribution of megaherbivores (S3 Fig). The Woodland cover and it's class metrics found to contribute less in explaining the distribution of megaherbivores and the percentage of deviation explained ranged from 1.9% to 2.8%.

Simulated future land cover for the year 2028

The future land cover for the year 2028 was simulated by adopting ANN algorithms, using the past temporal classified images of 2008 and 2018. The future land cover model was based on topographic functions. The percentage correctness and kappa (overall) was found to be 90.25% and 0.83, respectively. The comparison of model output with the classified image of 2018 predicts major change may happen in Woodland and Grassland cover types (Fig 5).

Furthermore, the model also showed a decline of about 6.23% from 2018 to 2028 in grassland. Whereas, an increment of about 6.09% has been predicted in woodland cover in 2028. Moreover due to no predicted significant change on other landcover types indicates that woodland cover may increase replacing grassland in the 2028.

Discussions

The GNP and its surrounding landscape is one of the last habitats of rhino in India and also home to two other major herbivores, i.e. elephant and gaur which is surrounded by a network of road, railway line and expanding tourism infrastructures [84,85]. As per the latest estimates, the rhino population in the landscape has increased over the years, from eight individuals in the early 1980s to 50 as per the Forest Department latest census (S4 Table). The increasing population of rhino and other large herbivores in the small landscape of the park is posing serious conservation and management issues, such as the increase in human-animal conflicts [19].

The analysis presented here demonstrates a significant change in land cover over the three decades, with a decrease in grassland and a concomitant increase in woodland (Fig 3 and S2 Table). The decline in grassland cover is a serious conservation and management threat to the mega-herbivores, particularly rhinos, which depend upon grasslands habitats for palatable grasses and forage [86–90]. Results from the predictor inspection showed that the grassland cover along with its class metrics remained the top predictors in explaining the spatial distribution of all three mega-herbivores viz., Rhino, Elephant and Gaur in GNP (Fig 4). The significant influence of grassland on megaherbivores of GNP can be attributed to several studies, suggesting substantial consumption of grasses in the diet of megaherbivores, especially for rhinos [21–25]. Hence, loss of grassland habitats will be challenging for the park authorities in maintaining the long term viability of the population of mega-herbivores. Further, it will greatly impact the overall carrying capacity of rhino, and other mega-herbivores in the GNP will decrease, as a result of the loss of grassland habitat [91]. Furthermore, rhinos, which are primarily grazers, may also experience increased competition with other large herbivores, many of which have broader dietary preferences [92].

Results shows that shrubland cover has declined from 10.53 km² from 1998 to 8.56 km² in 2008, followed by a slight increase from 2008 and 2018 (Table 2). Response from predictor inspection also suggests the contribution of shrubland and its metrics in shaping the spatial distribution, as AREA_AM_shrubland was about 5.6%, followed by euclidian distance to shrubland contributes about 4.3%, and both variables showed a strong positive association with distribution of megaherbivores. Thus, it can be opined that shrubland has provided a transitional area between the most suitable grasslands and relatively less suitable woodland for megaherbivores, especially for rhinos. Therefore megaherbivores of GNP can negotiate shrubland as movement corridors, and these shrub patches may compensate for the loss of grassland in the Gorumara up to some extent [93].

Landscape-level indices found to be useful in explaining the habitat quality, the present study has identified the related contribution of landscape-level variables, i.e. SHDI and SHEI as top landscape-level contributors with the positive response, which correlates with the reducing heterogeneity in the study landscape. So the reduction in the landscape diversity by more domination of woodland will eventually lead to less suitable habitat for megaherbivores. Apart from diversity indices ED, PD NP and AI are found to be in the top twenty variables in explaining the spatial distribution of megaherbivores. Influence of such variables has been supported by studies which have highlighted the significance of landscape configuration parameters on a diverse group of animals [50,51].

The patch aggregation, as well as interspersion and juxtaposition of different classes, indicate uniformity among the patches over the three decades. The resulted increasing uniformity in the patches can be attributed to plantation activities of the West Bengal State Gorumara Wildlife Division-II in the initial working plans. Further, the increase in NP and PD of woodland cover suggests homogenous with confined spatial distribution. The reduced values of the mean shape index from 1998 to 2018 of woodland patches indicate that woodland patches have become less complex over the three decades. Moreover, other indices such as IJI, SHEI, SHDI, ED and AI at class level in the study landscape indicate that during 2008 and 2018 much of the landscape have homogenized since similar types of patches are getting aggregated with reduced shape complexity. The decrease of NP and PD of grassland patches from 1998 to 2018 is a major concern for the conservation and management of mega-herbivores in GNP. We also found that among the anthropogenic predictors, rail infrastructure contributed the most and showed a negative relation in shaping the spatial distribution of megaherbivores in GNP, this could be related with the fact that in the recent past the intensity of megaherbivores mortality because of railways network has increased considerably [94].

Management implication for mega-herbivores

The change dynamics at both landscape and class level metrics in GNP indicate the transition of grassland to other land-use types. Which seems to be not suitable as GNP harbours a large population of sympatric megaherbivores, including elephant, rhino and gaur [25]. The prominent causes are monoculture plantations and ill-planned developmental activities, which mostly took place in the earlier working plans. Most of the conversion of natural habitat to Teak (Tectona grandis) and Jarul (Lagerstroemia speciosa) monoculture plantation took place in the areas during the VI and VII working plan before the establishment of the National Park. Such type of plantation contributes presently about 34% of the total plantation in GNP and associated forested areas [26]. The monoculture plantation of teaks does not favour the diversity of wild animals since it negatively impacts the undergrowth which is essential for the growth of fodder species [95].

Moreover, dense canopy allows very little sunlight to pass and thus inhibits the growth of light-demanding plants having browse value for ungulates. However, tree-dwellers like Giant Squirrels and certain birds may benefit from woodland habitat type, but may not support the species diversity in a broad spectrum and megaherbivores. Thus we suggest that old teak and other monoculture plantation should be kept in targeted areas on a rotation basis and canopy opening programs shall be restarted. These cleared areas should be considered for artificial regeneration of grassland to restore the carrying capacity of the ecosystem. However, the fodder plantation should not be single species oriented; instead, it should be done by using multiple species. In Jaldapara Wildlife sanctuary, a nearby protected area plantation of single grass species (Saccharum narenga) has impacted the overall faunal diversity composition and changing the plant community structure of the area [96]. We also recommend a systematic study on the feeding habitats of megaherbivores in GNP and surrounding landscape so that plantation preference should be given to plant species which are dominating the diet of megaherbivores of the landscape.

The relatively small area of GNP is supporting more than 900 gaurs and 50 rhino, despite a variety of anthropogenic disturbances and pressures (S3 Table). GNP has already exceeded the carrying capacity threshold for rhino, i.e. 43 individuals [19,37] and a further increase in the population will motivate animals to move out of the forested areas into the surrounding agricultural landscape in search of food [85,97]. Such movements are already causing human-wildlife conflicts in the landscape, and may antagonize the local human communities, which, in turn, may create a serious threat to the animals in question. One of the fundamental reasons for the uncontrolled increase in gaur population is food chain imbalance, in the near absence of top predators, i.e. tigers [85,97]. At present, there are no tigers left in the entire North Bengal landscape [98]. However, leopards are present in the landscape, but gaur is not their natural prey species [99]. This is leading to an ecological crisis in the landscape as rhinos are frequently found outside the park (far-flung as 50 km) to meet their foraging needs some of them staying seasonally or even permanently resulting in human-animal conflict in the landscape [15].

Furthermore, it will be imperative that tiger reintroduction in the landscape may be taken up by the forest managers for balancing the food chain. Alternatively, translocation of rhinos from GNP to other regions in the north Bengal will be a welcome step towards conservation and management of megaherbivores in the landscape. A similar initiative of rhino translocation from Pobitora Wildlife Sanctuary and Kaziranga National Park to Manas National Park through the Indian Rhino Vision 2020 (IRV 2020) program has accomplished the goal of conservation and population management of rhino [100]. We suggest establishing a rhino population in Buxa Tiger Reserve through shifting a few individuals from GNP. Apart from Buxa, the state government of West Bengal has identified a small reserve in Patlakhawa of Cooch Behar both located in the North West Bengal which can become a second home for GNP rhinos in the near future [101]. Moreover increasing connectivity by different habitat management program in the buffer forests of GNP along with increasing connectivity between Jaldapara-Gorumara and inclusion of territorial forests between Chapramari wildlife sanctuary (CWLS) and GNP in wildlife division II will be the most needed initiatives.

In addition, the forested area of GNP is crossed by an national highway (NH 31) as well as by railway lines, which is leading to loss of animals because of collisions with vehicles and trains [85,94]. These collisions of animals mostly occur on railways line between Siliguri and Alipurduar Junctions of about 150 km length. This railway line is known to traverse forest over 74 km about 44% of the length and passes through three protected areas including nine identified elephant sensitive movement corridors [94]. We recommend a reduction of train speed in this 74km stretch of the railway line to a maximum speed of 40km/hrs to avoid the collision, in conjunction with increased patrolling along the tracts jointly by railway and forest staffs for immediate implementation. Moreover, to adequately address this issue, construction of overpasses, underpasses, ecoducts and viaducts for providing safe passage at critical points and realignment of the present tract.

Limitations and methodological considerations

For projecting simulated future landcover for the year 2028, the cellular-Automata simulation was used based on Monte Carlo algorithm. The final simulated map for transition potential was based on classified images of 2008 and 2018. The models account only the previous change and not any anthropogenic or catastrophic natural processes. Minor false classification or mixed pixel classification in some instances may impact the classification and landscape heterogeneity [102]. The present study resolution of 30m from the satellite data (Landsat TM and OLT) may be considered low for certain landcover types.

Moreover, catastrophic effects like forest fires may also contribute to the landcover changes [103–106]. For evaluating the responses of different predictor variables for megaherbivores in GNP, we did not consider the distribution of other sympatric mesoherbivores, which might influence the spatial distribution of our study species. Because of the uncertainty introduced by both the landscape modelling and predictor response segments, authors emphasize that the present study represents an approach based on the available data and possible methodical considerations for megaherbivores in the study landscape.

Further efforts could improve our findings by including range and niche overlap between sympatric megaherbivores or by incorporating other plausible biotic/abiotic variable into the model environment. Despite these assumptions and caveats, the present study represents the best available postulate to evaluate the trend and effect of landcover change at a fine resolution. The trends of past and future prediction can be used to tackle a major management problem, from which the local authorities are struggling in recent times.

Supporting information

S1 Table [pdf]
Representing class-level and landscape-level metrics evaluated for landscape configuration change analysis for GNP.

S2 Table [a]
Accuracy assessment table for land cover classification of GNP for the year.

S3 Table [pdf]
Representing the transition in hectors from 1998–2008 and from 2008–2018.

S4 Table [pdf]
Representing the available population and demographic data for megaherbivores of Gorumara and North Bengal.

S1 Fig [pdf]
Landcover transition map of GNP.

S2 Fig [a]

S3 Fig [pdf]
Bar graph representing percentage of deviation explained by all covariates from Generalized Additive Modeling using the Software for Assisted Habitat Modeling (SAHM).

S4 Fig [pdf]
Methodological flow chart of the study.