Feasibility of reintroducing grassland megaherbivores, the greater one-horned rhinoceros, and swamp buffalo within their historic global range

Reintroduction of endangered species is an effective and increasingly important conservation strategy once threats have been addressed. The greater one-horned rhinoceros and swamp buffalo have declined through historic hunting and habitat loss. We identify and evaluate available habitat across their historic range (India, Nepal, and Bhutan) for reintroducing viable populations. We used Species Distribution Models in Maxent to identify potential habitats and evaluated model-identified sites through field visits, interviews of wildlife managers, literature, and population-habitat viability analysis. We prioritize sites based on size, quality, protection, management effectiveness, biotic pressures, and potential of conflict with communities. Our results suggest that populations greater than 50 for rhinoceros and 100 for buffalo were less susceptible to extinction, and could withstand some poaching, especially if supplemented or managed as a metapopulation. We note some reluctance by managers to reintroduce rhinoceros due to high costs associated with subsequent protection. Our analysis subsequently prioritised Corbett and Valmiki, for rhino reintroduction and transboundary complexes of Chitwan-Parsa-Valmiki and Dudhwa-Pilibhit-Shuklaphanta-Bardia for buffalo reintroductions. Establishing new safety-nets and supplementing existing populations of these megaherbivores would ensure their continued survival and harness their beneficial effect on ecosystems and conspecifics like pygmy hog, hispid hare, swamp deer, hog deer, and Bengal florican.

Species are facing an unprecedented extinction crisis due to anthropogenic impacts 1 with large carnivores and megaherbivores bearing the brunt 2 .These taxa have been extirpated from many of their natural habitats by direct hunting for meat, trophies, crop protection, and retaliatory killing 3,4 . However, large carnivores have been observed to recolonize areas where these threats have been removed, if such areas are connected with source populations 5 . Habitat connectivity for megaherbivores differs to connectivity in carnivores, since carnivores can often pass through degraded habitats, while megaherbivore corridors require forage, water and absence form potential conflict with humans. Megaherbivores share slow life history traits, which make them more vulnerable to threats of habitat loss and poaching 6 . As a consequence, recolonization is rarely observed in megaherbivores other than in elephants 7 , since most source populations are depleted and connecting corridors are degraded or destroyed. Planned reintroductions can be a useful strategy to establish safety-net populations and are often a better option, since natural recolonization is a slow and stochastic process 8 .
Megaherbivore species are largely unaffected by non-human predation, so can sustain high stable densities, allowing them to greatly modify the ecosystems they inhabit 9 . As ecosystem engineers, megaherbivores regulate tall coarse grass through grazing and trampling, which allows the growth of palatable species accessible for consumption by smaller mesoherbivore species 10,11 . Reintroduction and supplementation of megaherbivores would not only create safety-net populations but also restore an important ecological role in currently degraded habitats 12 and recover the potential of such habitats to sustain historical faunal assemblages.

Results
Species distribution modeling. Greater one-horned rhinoceros. Based on Area Under the Curve (AUC) of the Receiver Operator Curve (ROC) and on 100 bootstrap runs of omission/commission analysis with 20% test data showed that the rhinoceros model had a good fit [AUC = 0.96 (SE 0.0007)] and predictive ability (Fig. 1a,b). This model was also rated the best by True Skill Statistics (TSS) and had the lowest Akaike Information Criteria (AIC). Rhinoceros occurrence was best explained by (a) Distance from grassland, (b) Distance from forest, (c) Maximum temperature of warmest month (d) Annual Precipitation and, (e) Distance from water ( Table 1). As expected, species response curve showed a decline in habitat suitability with increasing distance from grasslands (Fig. 1c). Rhinoceros were unlikely to occur at distances > 2 km from grasslands. Distance from grassland habitat explained the maximum variation (69%) in the rhinoceros distribution model. Rhinoceros occurrence declines rapidly from forest edge (Fig. 1d) which contributed 16% to the model. A temperature range between 32 to 40 °C during the warmest months of the year governed occurrence of rhinoceros (Fig. 1e). This covariate (Maximum temperature of the warmest month) explained 8.2% of the variation in the data in the best model. Rhinoceros www.nature.com/scientificreports/ occurred where annual rainfall exceeded 1700 mm and the probability of occurrence increased with increase in rainfall up to 4000 mm (Fig. 1f). Rhinoceros occurrence probability declined sharply with increasing distance to water up to 5 km, after this distance there was high variability in the response curve, This covariate contributed the least to the model at 1.9% (Fig. 1g).
Swamp buffalo. The Buffalo model had a good fit with an AUC of 0.95 (SE 0.0027) but the predictive ability based on 100 bootstrap runs of omission/commission analysis with 20% test data (Fig. 2a,b) was not as good as that for rhinoceros. The TSS was the best for this model while the AIC was second best (Table 1). However, since this model had a good fit and made ecological sense supported by TSS we use it as the best model. The covariates of the best buffalo model were (a) Distance from grassland, (b) Distance from forest, (c) Annual Precipitation and, d) Distance from water (Table 1). Distance to grassland had the highest contribution to the model. Buffalo occurred in the proximity of grasslands and were not likely to be found beyond 1 km distance from grasslands (Fig. 2c). The model showed high buffalo habitat suitability at forest edges with low probability with increasing distances from forests (Fig. 2d). Areas having annual rainfall above 1500 mm were preferred (Fig. 2e). Habitat suitability for buffalo declined rapidly with increasing distance from water for up to 2 km after which habitat was unsuitable with high variability in model predictions (Fig. 2f).  Table S1). Medium sized populations (K = 20-30) are shown to be viable when initial reintroductions are undertaken with > 8-10 rhinoceros and sites are occasionally supplemented, but these populations would not withstand poaching ( Table 2, Scenarios 5-11 and Table S1,). Populations with K ≥ 50 have better chances of survival which increases with supplementation. These populations would also withstand low levels of poaching when supplemented at initial stages ( Table 2, Scenarios 14 & 15). Populations ~ 100 would survive long-term, even when subjected to a low-level poaching, and are able to retain high levels of heterozygosity even without immigrants ( Table 2, Scenario 19-21  and Table S1). Populations with K ≥ 100 are ideal for long-term persistence ( Table 2, Scenario 22). A metapopulation comprising of Kaziranga, Orang and Laokhowa-Bura Chapori demonstrates a stochastic growth rate of 0.022 (SE 0.019). The extinction probability of that metapopulation is zero and heterozygosity is maintained at ~ 100%.
Swamp buffalo. Small populations (K ≤ 20) exhibit low persistence probability despite supplementation ( Table 3, Scenarios 1-4 and TableS2). Medium sized populations (K = 50) could persist with initial reintroduction of > 10 individuals supplemented for a decade, but would not sustain poaching offtake (Table 3, Scenarios 5-9 and Table S2). Large populations (K = 100-200 buffalos) remain viable in settings which experience natural catastrophes, and are also resilient to moderate poaching losses if initial founding population is > 30. Populations are sensitive to the size of founding population and depend upon continued supplementation (Table 3, Scenario 10-16 and Table S2). Areas which could sustain (K) > 250 were ideal for long-term persistence and could tolerate moderate poaching. However, all populations were sensitive to founding population size, and a founding population > 30 is optimal.
Identifying and prioritizing suitable habitats. The species distribution probability asc. layer obtained as the median from 100 Maxent bootstrap runs were exported to Arc GIS 10.5 to produce probability maps for rhinoceros ( Fig. 3) and buffalo (Fig. 4), showed reasonable extents of suitable habitat outside of these species' current range. Maps from conservative estimates of 95% lower limits also identified substantial patches of suitable habitat for reintroductions ( Fig S1)

Discussion
Reintroductions offer great potential for species recovery when following scientifically-informed strategies and relevant risk assessment 22 . Our evaluation of potential reintroduction sites and extant populations addresses these aspects, although any subsequent site-specific plans must also consider IUCN guidelines. Several identified patches sit inside Protected Areas (PAs) with some ready to receive animals, since threats have been addressed and park management capacity exists. We highlight the most feasible actions needed at potential sites. Reintroduction would create safety-net populations for megaherbivores and restore their ecosystem-engineering role, benefitting endangered grassland specialists including pygmy hog (Porcula salvania), hispid hare (Caprolagus hispidus), swamp deer (Rucervus duvaucelii), hog deer (Axis porcinus), and Bengal florican (Houbaropsis bengalensis).
We used a well-established modeling approach of Maxent models 23 with presence locations obtained from all extant populations to identify potential habitats. Since the modelling space defined in Maxent was clipped within limits of current species' distribution (elevation, temperature, precipitation), the response curves of these variables should be viewed accordingly; we did not model suitable range in Maxent, but rather the high preference habitats within this suitable range. Additionally, covariates used in the final model were uncorrelated, so the response curves of each covariate are similar when compared to individual and cumulative contribution with all other covariates (Fig. S2). By capturing variability in model-based inference on 100 bootstrap runs for 95% lower limits, inferences remain conservative.
All our presence locations were within PAs (since surviving populations exist only in PAs), so the metric 'distance to PA' would overpredict potential for occurrence in most PAs, we therefore refrained from using this covariate in our final models. Instead we accounted for legal status of habitats after the SDM modelling process during decision-making. Human Footprint Index did not conform to a priori hypothesis and was not used Table 2. Results of 500 simulations of population trajectories over 100 years in VORTEX (9.93) to assess the viability of greater one-horned rhinoceros populations with different scenarios of carrying capacity, poaching, catastrophe, initial population size and supplementation. Here AF -adult female, AM -adult male, r = growth rate of population, (SD) = standard deviation, N = population size at the end of 100 years, PE = probability of Extinction and H = heterozygosity of the population at the end of 100 years.

Scenario
Carrying capacity Initial population Supplementation Frequency of catastrophes Frequency of harvest r (SD) PE N H% www.nature.com/scientificreports/ after exploratory analysis, since most of the fertile floodplain habitats outside the realm of legal protection were already occupied for human land use. Species presence locations were often juxtaposed with pixels having high Human Footprint values due to hard PA boundaries, making it ecologically uninformative. Variables such as 'distance from natural grasslands' and 'distance from forest' better captured negative influences of human impact compared to Human Footprint Index and were aligned with a priori expectations.
Often the test data of species occurrence used to validate models are in close proximity to species presence points used for model building, causing spatial autocorrelation 24 that inflates the accuracy (AUC) of the model 25 .
In both megaherbivore models geographic clustering of training and test data may have caused spatial autocorrelation and inflation of AUC values. However, our assertion is that species occurrence data used in our models spans the entire extant geographic range of both megaherbivores across the entire spectrum of covariate space occupied by each species. We use only one presence location per one km 2 pixel (as explained by Philips et al. 2017 85 ), and use the bias correction file option in Maxent to address this limitation of training data 26 . We also use TSS and AIC to evaluate our models.
Our habitat suitability analysis showed PAs such as Kishanpur WLS, Sohagibarwa WLS and Nameri TR had suitable habitat to support ≤ 10 rhinoceros. Results of PHVA for rhinoceros suggest small populations (≤ 10) had low probability of persistence so these areas are rejected. Suitable habitat identified in Rajaji NP was in the Table 3. Results of 500 runs of population habitat viability analysis of wild swamp buffalo populations with different scenarios of carrying capacity, poaching, catastrophe, initial population size and supplementation over a period of 100 years. Here AF -adult female, AM-adult male, r = growth rate of population, (SD) = standard deviation, N = population size at the end of 100 years, PE = probability of Extinction and H = heterozygosity of the population at the end of 100 years.

Scenario
Carrying capacity Initial population Supplementation www.nature.com/scientificreports/ floodplains of the Ganges and its tributaries, an area heavily utilized for religious pilgrimage in proximity to townships of Haridwar and Rishikesh, again unsuitable for reintroduction due to potential human-wildlife conflict. Although our analysis suggests that Valmiki itself could support only a small population (35-40 individuals), we recommend reintroduction, since Valmiki receives natural immigrants from Chitwan (Nepal), so offers probability of long-term metapopulation persistence. That said, investment to re-route the railway line passing through prime rhinoceros habitat is necessary, and is being considered by the Bihar Government (Valmiki Park Director 2018, Personal Communication). An increase in law enforcement and reduction of human activities are also required, although Valmiki has clearly improved its protection regime in recent years as evidenced by an increasing tiger population 27 , a species vulnerable to poaching 28 .
Encouragingly, existing rhino ranges in Shuklaphanta, Bardia, and Orang have potential to sustain higher numbers. Supplementation would be possible if there is investment in protection, habitat management, and reduction of anthropogenic pressures and human-rhinoceros conflict. The same is needed for existing     (Table 4). Site evaluation and habitat studies (including grassland food-source species like Saccrum, Imparata, Arundo and Vetivaria) 31 suggest that Corbett is suitable, with ample surface water (rivers, reservoirs, pools) 32 plus sufficient grassland patches interspersed with forest cover, to support mixed foraging and cover during calving and in winters 33 . Habitat management such as artificial wallows may be required, since the largely Bhabhar landscape lacks muddy swampland and oxbow lakes of the Terai. The level of protection in the park is high with managers capable of handling reintroduction 34 . A key benefit of Corbett over other PAs is political stability in Uttarakhand, including in and around the PA (in contrast to Assam and Nepal 18 ). The large valleys of Corbett TR offer ideal habitats for rhinoceros, bounded by the lower Himalayas (north), Shivalik Hills (south) and the Ramganga Reservoir, retaining rhinoceros within the PA, so limiting conflict.
Laukhowa-Bura-Chapori complex, proposed for reintroductions by Rhino Vision 2020 36 , was also identified by our study, however considerable conservation investments are needed before reintroduction (Fig. S3). Major interventions needed include (a) reduction of cattle grazing, (b) grassland management and weed control, (c) protection (guards, weapons, law enforcement training, vehicles, and infrastructure). The site is of particularly significance 27 , as it connects habitat along the Brahmaputra river for rhinoceros, elephants, buffalo and tigers dispersing from Kaziranga. Our PHVA analysis shows that if rhinoceros in Orang TR, Kaziranga NP and Laokhowa-Bura Chapori WLS complex are managed as a metapopulation they will persist even in the small reserves of Orang and Laokhowa-Bura Chapori for the long-term. Although Dibru Saikhowa WLS (Rhino Vision 2020 site) 36 and D'Ering offer extensive habitat, both are affected by seasonal agriculture encroachment and high livestock density. Effort to remove or reduce those threats would enable these island sanctuaries to hold an estimated > 100 rhinoceros each.
One potential metapopulation includes the complex of Bardia and Shuklaphanta of Nepal across the international border to the Indian complex of Dudhwa, Katerniagath, and Pilibhit along the Khata Corridor, Lagga-Bagga WLS and flood-plains of Sharda river 37 . Rhinoceros sometimes transit these PAs, however new proposed border roads will create further obstacles to transboundary wildlife movements. Metapopulation connectivity would benefit other species like (tigers, elephants, buffalos), reducing the need for supplementation.
Corbett is isolated with no adjacent rhinoceros or buffalo populations, however its large carrying capacity would make a population viable if genetically diverse founders are selected, supported by supplementation. Hastinapur WLS in Uttar Pradesh 38 and Surai Range in Uttarakhand were previously considered for rhinoceros, but our analysis rejects both. Hastinapur WLS lacks forest cover (essential for thermoregulation and calving) and is bordered by agriculture and high human density. The Surai range, although Terai habitat, is small and suffers human encroachment.
We have modelled supplementation for a period of 5-10 years, a time period coinciding with most government funded projects. If however, supplementation continues intermittently over the long-term it would dramatically improve persistence and genetic variability of populations. Our population estimates from model-based computations for Chitwan and Kaziranga match the current estimates of rhinoceros at these sites 15 which suggests that extant populations are nearing carrying capacity and can serve as source populations 39 . Pobitora WLS, which has close to 100 rhinoceros and a sizable buffalo population, was not identified as good habitat in our results. This shows the limitations of model-driven inferences, which may fail to account for outliers resulting from species' behavioral plasticity or exceptional tolerance by local people. Pobitora is unique, in being imbedded in a human dominated landscape with minimal availability of optimal conditions, yet rhinoceros and buffalo persist in this remnant (39 km 2 ) due to community tolerance (Fig. S5).
Our results show (somewhat unexpectedly considering illegal demand for rhinoceros horn) that buffalo are more prone to extinction than rhinoceros, which is actually congruent with the IUCN Red List 40 which categorises the rhinoceros as vulnerable and buffalo as endangered. Buffalo are rarely poached commercially, but our results suggest that even small off-take due to bush-meat demands would be a major cause of concern. Our PHVA results for buffalo are likely conservative since studies on demographic parameters of wild buffalo are sparse and potentially biased towards higher mortality estimates.
Many buffalo populations of Northeast India are either too small or possibly hybridized with domestic buffalo (derived from swamp buffalo) 16 . In western regions buffalo are largely locally extinct including in locations which could potentially sustain populations. It would be prudent to reintroduce Bubalus arnee into suitable western sites, since domestic buffalo in those localities are derived from river buffalos (Bubalus bubalus) 41,42 , presenting less risk of hybridization.
Priority sites for buffalo reintroduction in the subcontinent are Chitwan NP-Valmiki TR complex and Shuklaphanta-Bardia-Dudhwa-Katerniagath-Pilibhit complex since both could sustain > 400 individuals. PAs like Sonai Rupai, D'Ering, Dibru Saikhowa, Jaldapara, and Royal Manas could each hold more than 100 individuals, with the potential to be managed as metapopulations. However, bushmeat hunting persists in northeastern India 43 , which poses a major threat to wild buffalo. Higher densities of buffalo may not be possible until human pressures are reduced and management of invasive plants is undertaken. The model-inferred median population of buffalo sustainable in Corbett NP is about 130 individuals. However, there was high level of uncertainty in our model-based estimates for Corbett (Table 4), which we believe was due to the sampling of presence locations being limited to extant populations in the far east region, which experiences different bioclimatic conditions to Corbett in the western Terai. www.nature.com/scientificreports/ The reintroduction of buffalo could act as a test-bed for rhinoceros reintroductions, enabling establishment of management protocols, enhancing staff experience and ensuring effective future operations for rhinoceros. Reintroduction of buffalo to Dudhwa should be feasible since park managers already have protection and monitoring systems in place following experience of managing reintroduced rhinoceros. Chitwan NP in Chitwan-Parsa Complex of Nepal (a high potential site for buffalo reintroduction in this study) has seen recent buffalo reintroduction, which has a good chance of establishment due to habitat quality, local people's awareness of ecosystem services, experience of megaherbivores and economic benefits from wildlife tourism 44 .
Knowledge of population genetics must inform any reintroduction. For buffalos, hybridization with domestic buffalo is a key factor, and animals of pure wild lineage are needed as founders for reintroduced populations. Possibly the last pure representatives of Bubalus arnee are around 50 buffalo surviving in central India (Udanti-Sitanadi TR, Indravati, in Chhattisgarh and Sironcha in Maharashtra) 45,46 . According to our PHVA models, this Central Indian population suffers a high probability of extinction. The animals reside in a zone of high political unrest, making local conservation difficult 47 . Efforts to rescue these buffalo is vital, including taking animals into captivity for conservation breeding, or translocation.
An important aspect highlighted by our SDM was that both the greater one-horned rhinoceros and the wild swamp buffalo have very narrow tolerance of temperature and precipitation (Fig. 1e,f; S5). This has serious implications for persistence of populations under climate change. Current populations of both species are limited to PAs largely surrounded by human land-use, making natural range-shifts unlikely 48 , so conservation policy and strategy must address this limitation by managing metapopulations. Western (Corbett NP) and central terai (Dudhwa-Katerniagath-Bardia) are at higher elevations compared to the Brahmaputra plains and could provide refuges in the advent of climate change.
Several factors limit the potential for reintroducing species in their historic range, including loss and degradation of habitat and competition with human interests. For rhinoceros, even in areas with good habitat and least potential for conflict with humans, managers were reluctant to reintroduce the species, due to resource demands and likely media and political pressure in the face of potential losses of animals. Suitable habitats and management are insufficient for protecting species with high value in illegal wildlife markets; only changes in human perceptions, values, and strict law enforcement will address illegal demand 49 .
We advocate the establishment of safety-net populations across a varied geographical range, and our recommendations address the uncertainties associated with climate change. Continued anthropogenic pressures mean that maintaining effective ecosystems will require balanced species assemblages. Reintroduction and supplementation of viable populations of keystone species are important tools for sustaining and enhancing vital functions in critical natural systems.

Materials and method
Field sampling. Field sampling was undertaken to obtain extant rhinoceros and buffalo locations (from India, Nepal and Bhutan) for modelling habitat suitability and assess the reintroduction potential of identified sites. The current extant and potential reintroduction sites were assessed for level of protection, anthropogenic pressures, size, and quality of available habitat for the target species through site visits and interviews with wildlife managers. The interviews addressed managers' perception of management strengths and weaknesses, and attitude towards reintroduction/supplementation of both megaherbivores (Table S3). All interviews were conducted after obtaining informed consent of the respondents. Prior ethical approval for the study was granted by the School of Anthropology and Conservation, University of Kent, UK, (Reference ID-48-PGT-17/18) following the standards of the American Anthropological Association. All research was conducted under the laws, rules, and regulations of each country. Interviews were conducted in English since all wildlife managers interviewed were conversant with English. Information from site visits and interviews were used to evaluate sites that were selected by SDMs for practicality of reintroductions and potential establishment of populations.
Species distribution modelling. Several approaches to modelling SDM are available for data typ used in this study i.e. Presence vs. Background points. These include like (1) Ecological niche factor analysis (ENFA) 50 , (2) Genetic algorithm for Rule-set prediction (GARP) 51 (3) Maximum Entropy models (Maxent) 53 . Earlier researchers have compared the performance of these approaches and found that Maxent models outperform the other approaches 52 . We used Maximum Entropy Species Distribution Modelling in Maxent (Version 3.4.1) 53,54 for modelling potential habitat of rhinoceros and buffalo . Maxent uses machine learning to develop relationships from known species occurrence and background data with ecologically meaningful spatial environmental covariates 54 . The program subsequently predicts potential distribution of species from covariate relationships across modelled space 53 .
Species occurrence data. Species occurrence information from extant populations of both species was obtained by direct sightings with coordinates recorded by hand-held GPS, camera trap photo-capture locations across the Terai-Brahmaputra floodplains in India 27 , and from locations of satellite-GPS-radio-collared rhinoceros in Chitwan NP 14 . To avoid spatial autocorrelation and oversampling, we picked only one location from an area of 1 km 255 . A total of 358 rhinoceros locations and 78 buffalo locations were used for Maxent modelling. Historical records of occurrence were not used since many covariates used in our models have changed since historical times. Furthermore, our modelling objective was to identify extant available habitats and not historical species range.
Eco-geographical variables. The following covariates were used as predictive variables for SDM: Distance from grassland www.nature.com/scientificreports/ Grasslands of the study area were digitized using unsupervised and supervised classification 56 in ESRI ArcMap 10.5.1 (Environmental System Research Institute 2016) using LandSat 8 imagery 57 and corrected using Google Earth. We ground-validated grassland locations by obtaining coordinates from a hand-held GPS device (Fig S6). Of the 141 ground validation points 125 were correctly classified resulting in an accuracy of 88.6%. Euclidean Distance tool in ArcMap 10.5.1 was used to calculate the distance of each pixel from grassland polygons at a resolution of 30 m. Both rhinoceros and buffalo can be considered as grassland dependent species 46,58 . A priori we expected both rhinoceros and buffalo occurrence probability to be high within, and on the edge of, grassland habitats and to decline rapidly with increasing distance to grasslands.
Distance from forest Data on forest cover was obtained from GlobCover (2009) at a spatial resolution of 300m 59 . Euclidean Distance of each pixel to a forest patch was computed using ESRI ArcMap (10.5.1). Both study species inhabit grassland-forest mosaic habitats 46,58 . We expected both species to have high occurrence in forest-edge habitat with a declining response for increased distance to forests.
Distance from water Spectral indicators such as Normalized Difference Water Index and Modified Normalized Difference Water index were used to classify waterbodies (> 30m 2 ) from LandSat 8 imagery 60 . However, smaller water bodies which are also used by these species could not be detected by satellite imagery and therefore not accounted for in our models. Euclidean Distance to water for each pixel was calculated using ArcMap 10.5.1. Both rhinoceros and buffalo need ample access to surface water and spend a lot of their time wallowing for thermoregulation 61 . In Subedi's 33 study all radio-collared rhinoceros' locations in Chitwan NP were within 1.8 km distance from water bodies. An essential element of rhinoceros and buffalo habitats are flowing rivers, pools and oxbow lakes. Studies on behaviour of both species show the use of pools and oxbow lakes for wallowing especially in the dry season 45,62 . Our expectation was high occurrence in proximity to water with a rapid decline in species occurrence with increasing distance from surface water.
Normalised Difference Vegetation Index (NDVI) NDVI composites were derived from Moderate Resolution Imaging Spectroradiometer (MODIS) data, which was obtained from online data pool, developed by NASA Land Processes Distributed Active Archive Centres. MODIS products are available at spatial resolution of 250 m with 16-day interval cycle. Use of NDVI for species modelling enable better predictions when used alongside other ecological variables 63 . In this study, three values of NDVI were used (1) Pre-Monsoon NDVI (March-April), which indexes canopy cover of the driest months, (2) Post-Monsoon NDVI (October-November) which characterise maximum cover, and (3) Difference in NDVI of Post-Monsoon and Pre-Monsoon which indicates monsoonal flush of annuals compared to more perennial canopy cover. We expected rhinoceros and buffalo to have a positive response in occurrence probability to higher pre-monsoon NDVI as these pixels would be correlated with forage availability in resource lean summers. Species occurrence response to NDVI difference was expected to increase with increasing difference, as higher difference to reflects fresh flush of nutritive forage 64,65 .
Geo-climatic data Climatic data are commonly used for modelling species distribution 60,61,66 . Both species prefer wetter habitats and moderate climate 33,46 . The modelling extent was clipped by the limits of species occurrence recorded for elevation, temperature and rainfall 67 . Data on (i) annual rainfall (b15, layer clipped for annual rainfall above 1000 mm), (ii) precipitation of the driest quarter (b17), (iii) maximum temperature of the hottest month (b4, layer clipped below 45 °C temperature), (iv) minimum temperature of the coldest month (b6, layer clipped above 0 °C temperatures) and (v) precipitation of the wettest month (b13) were obtained from Worldclim website 68 at a resolution of 0.5 km 2 . We hypothesized that both species would show a parabolic or peaked response to rainfall and temperature, with a skew toward higher rainfall and a peak at moderate temperatures.
Elevation plays a major role in the distribution of both the megaherbivores as they prefer flatter and less rugged terrain 9 . Ecological studies on rhinoceros indicate preferences for low-lying flat habitats 58 . However, for buffalo, studies have recorded presence up to 1000 m elevation, so we clipped our modelling extent to elevation below 1000 m. We expected a sharp decline in rhinoceros occurrence with increasing elevation and a similar but shallower decline for buffalo. We used the Digital-Elevation-Model produced by the Shuttle Radar Topography Mission (a joint endeavour by NASA, National Geospatial-intelligence Agency, and German and Italian Agencies) 69 .
Distance from protected area (PA) Both megaherbivores suffer poaching and their current populations are largely restricted to PAs 35,46 . We expected species occurrences to decline as distance from PA increased. PA boundaries were taken from the Protected Planet website (https ://www.prote ctedp lanet .net/), the Wildlife database cell of the Wildlife Institute of India and the Project Tiger Directorate. The distance of each pixel from PA was calculated using the Euclidean Distance Tool in ESRI ArcMap (10.5.1).
Human footprint Megaherbivores often conflict with humans due their requirement of vast fertile plains, propensity to crop raiding, damage to human property and lives 9 . Proximity to human settlements also raises the problems of habitat degradation, encroachment, poaching and hunting 70 . We therefore, expected a negative relationship between human footprint index and megaherbivore presence. The Global Human Footprint Index dataset was obtained from Last of the Wild Project 71 , which is the Human Influence Index normalized by biome and realm.
Modelling in maxent. The habitat covariates were exported to ArcMap (10.2) to obtain uniform resolution and concordance with grids of 0. 83  www.nature.com/scientificreports/ computed in Eris ArcMap (10.5.1) ( Table S4). The correlation threshold was kept at r > ± 0.70 to check for redundancy in information. No two correlated variables were used together in a single model. Maxent operates by developing relationships between known occurrence locations and covariates by comparing them with background locations. Species locations used for training our models were from a restricted part of the modelled space (though all extant species locations were sampled). We therefore used a bias correction file 26 created using inverse distance weight (IDW) interpolation method in ArcMap (10.5.1) to guide Maxent into picking background locations from space containing occurrence locations 54,72 . This approach would enhance the accuracy of our models by limiting the predictive relationships to be developed from the extent of presence vs background from the same area. We used 80% of presence locations of each species to train Maxent models and the remaining 20% were used to test models. The run type was set at Auto with Linear, Hinge, Quadratic and Product functions were selected in Maxent to model relationships between covariates and species occurrence. The maximum number of background points was set at 10,000. Regularization Multiplier was set at default (1) for both species models. We first explored univariate relationships between covariates and species occurrence. Subsequently, based on the univariate responses, we used the covariate combinations that were ecologically meaningful and represented climate, abiotic and biotic habitat and human disturbance to develop more complex models for our species distribution models.
One hundred bootstrap simulations were run for the best model for both species. The median prediction was used to evaluate sites while 95% upper and lower predictions were used to capture the variability associated with model uncertainty. We report species occurrence probability as spatial maps. To compute patch-areas suitable for reintroduction, we used 'average of maximum training sensitivity plus specificity cumulative threshold' to categorise pixels as suitable habitats 23 . We computed 95% confidence intervals on suitable habitats using the threshold value. Covariate selection and model evaluation was based on (a) Receiver Operation Characteristic (ROC) curve 73 , (b) contribution of covariates individually and in ecologically meaningful combinations to explain the variability of the training and test data sets and (c) Omission/Commission analysis of test dataset (d) True Skill Statistics (TSS) 74,75 and (e) Akaike Information Criteria (AIC) 52,76 Species response curves to each covariate were examined and ecologically interpreted 77,78 .
Rhinoceros and buffalo were observed to be at high densities in Kaziranga (6 rhinoceros 79 /km 2 ; and 3.2 buffalos/km 2 , https ://www.kazir anga-natio nal-park.com/wild-buffa loes.shtml ) while Chitwan had lower densities of rhinoceros (1.4 Rhinoceros/km 2 ) 80 . Gorumara NP in West Bengal includes forest, plantation and patchily distributed grasslands carrying a density of about 0.65 rhinoceros in the PA 81,82 . These differences were likely due to difference in productivity and habitat availability in these PAs; with Kaziranga being more productive and the entire PA being favourable habitat while Chitwan (constituted by Terai, Bhabhar, and Churia hills) and Gorumara having only parts of their habitat as favourable. We report crude densities for these PA but use a conservative estimate of one individual km -2 (and 2 km -2 for Orang TR which has habitat similar to Kaziranga) as the potential ecological density of both species for reintroduction sites. Achievable population sizes were computed by multiplying only the suitable habitat area by this density.

Population habitat viability analysis (PHVA).
We assessed the viability of potential reintroduced and extant populations using PHVA in Vortex 9.93 21 . PHVA allows assessment of persistence of a population over a specified timespan incorporating the stochastic nature of demographic parameters, carrying capacity, environmental fluctuations, genetic variability, and deterministic impacts such as poaching 83 . Information on demographic parameters for rhinoceros was obtained from previous studies in Chitwan 39,58,67 . Demographic parameters for the buffalo were also obtained from literature 16,84 (Table S5). We modelled populations with carrying capacities (K) ranging between 10-150 for rhinoceros and 20-500 for buffalos to reflect habitat patch sizes estimated through Maxent. Each population was subjected to scenarios with varying sizes of initial population, levels of supplementation, poaching and mortality caused by catastrophe (floods/diseases). Based on population sizes and potential movement between Kaziranga, Orang and Laukhowa-Burachapori, a metapopulation of rhinoceros was modelled. This model provided information on the role of movement across the landscape upon viability of individual populations and the metapopulation. Five hundred simulations for each scenario were run for 100 years. Extinction probability, stochastic growth rate, population size, persistence time and the level of heterozygosity were used as evaluation criteria 39 for prioritizing potential sites.