Wintering habitat modelling for conservation of Eurasian vultures in northern India

Eurasian Black Vulture (EBV) and Eurasian Griffon Vulture (EGV), while residents elsewhere, winter in Uttar Pradesh, India. Knowledge of the habitat and regulating factors is obligatory for protection and better management of these vultures. Therefore, different types of habitats were mapped using eight species distribution models. Presence records from field survey, published data and citizen science, and 23 bioenvironmental raster layers were the model inputs. Eighteen models were developed whose strength varied greatly. As per the performance indicators, GBM and GLM were found to be superior models for EGV. For EBV all models were acceptable. MARS, with good model strength, was rejected on the grounds of field verification. However, the Ensemble model, overall, was found the best. As per this model, good habitat was restricted mostly in the Tarai ecozone. The top two vital variables were NDVI, and bio13 for both the vultures. The most vital temperature variable for EGV was bio08 while bio09 for EBV. Tarai ecozone showed the largest expanse of suitable area for both the vultures followed by Vindhyan-Bundelkhand, Gangetic plains and Semiarid ecozones. Among the two, EBV (49000 km) had more suitable area than EGV (37000 km). Agricultural areas were found to be largely unsuitable. As per land cover, good habitat was mostly confined in forests. For better management of these wintering vultures which need only roosting and foraging, it is proposed that destruction of forested habitat and decrease in foraging materials needed immediate attention and control.


Introduction Introduction Introduction Introduction
Vultures are well recognised scavengers providing vital ecosystem services of keeping the environment clean. Any threat to them leads to a disbalance of ecological equilibrium carrying risks of pollution and diseases in wildlife, humans and livestock (Jha and Jha, 2021b). This service is further boosted in the wintering habitat by the arrival of migratory vultures. Two of the three migratory vultures being harboured in India are of Eurasian origin. These two Eurasian vultures, Eurasian Black Vulture or Cinereous Vulture (Aegypius monachus Linnaeus 1766) hereafter, Eurasian Black Vulture (EBV) and Eurasian Griffon Vulture or Griffon Vulture (Gyps fulvus Hablizl 1783), hereafter Eurasian Griffon Vulture (EGV) (Figure 1) are common sight in north-west-central India during winters (Jha and Jha 2021a). Eurasian Black Vulture has a "near threatened" (NT) conservation status globally with a decreasing population trend. This species breeds in many countries of Europe and Asia and their wintering range extends to several other Asian countries including India (BirdLife International. 2021a). On the other hand, Eurasian Griffon Vulture has a "least concern" (LC) conservation status with an increasing population trend and also has a large breeding range in European and Asian countries. This species too migrates to spend the non-breeding season in northern African and Eurasian countries including India (BirdLife International, 2021b). Northern India is a well recognised wintering ground for these two vultures (Jha, 2015;MoEFCC, 2020). Wintering in simple terms is spending a good amount of time in better living conditions as compared to the original home of a bird. Bird migration is the regular seasonal movement between breeding and wintering grounds. Birds also show fidelity to wintering sites and visit them year after year (Fox et al., 1994;Clausen et al., 2018) until altered severely. This is one of the strategies of the birds to avoid harsh winter at its breeding grounds and survive for long. This indicates that individual sites have to be beneficial in terms of food abundance, sheltering structure, minimised disturbances and optimum weather conditions to survive (Anderson et al., 1992;Robertson and Cooke, 1999;Clausen et al., 2018). However, threats to migratory birds have grown with habitat destruction especially at stopovers and wintering sites (Yong et al., 2021). Therefore, it is pertinent to know about the available habitat expanse and conditions in order to strategize their conservation.
Delimiting habitat expanse, assessing suitable habitat and identifying the habitat variables are important components of conservation ecology which help strategize the action plan for recovery, maintenance and natural augmentation of avian population (Banda and Tassie, 2018;Nursamsi et al., 2018;D'Addario et al., 2019;Habibzadeh and Ludwig 2019;Kanth et al., 2020). Ecological niche or species distribution models provide immense opportunities to identify suitable habitats and bioenvironmental factors for current as well as future scenarios (Osborne and Seddon, 2012;D'Elia et al., 2015). Information generated by these models can be utilised for conservation planning and management purposes (Rodriguez et al., 2007;Guisan et al., 2013).
Coincidentally, there is a growing interest in utilizing habitat suitability models in spatial conservation planning for vultures in different parts of the world (D'Elia et al., 2015;Phipps et al., 2017;Farashi and Alizadeh-Noughani, 2018;Bahadur et al., 2019;Jha and Jha 2021a). Vector Machine (SVM), etc. (Li and Wang, 2013;Bucklin et al., 2015;Aguillar et al., 2016;Kafas et al., 2016;Jha, 2021). These habitat prediction algorithms of different orders like regression, classification technique and machine learning/complex models (Li and Wang, 2013), vary in prediction results for a particular species and environment (Konowalik and Nosol, 2021). No modelling technique is best for all situations (Marmion et al., 2009). Researchers proposed that rather than using single modelling technique, different models for each species should be performed and the most accurate technique should be selected (Heikkinen et al., 2006;Raina et al., 2014). However, Dormann (2008) concluded that selection of the single-best model introduces model selection error, which can be, to some extent, overcome by the multi-model approach. Several ecologists in general have adopted the strategy to develop an ensemble model of various best performing tools as it's components since it addresses the weakness of each individual algorithm and accommodates the positives of all the component algorithms (Stohlgren et al., 2010;Aguillar et al., 2016;Chalghaf et al., 2018;Früh et al., 2018;Ahmed et al., 2021).
In light of the above, the following objectives were set for this study on wintering Eurasian vultures of northern India to: (i) develop various species distribution based models and find the most suitable ones, (ii) project habitat area suitability under different categories, (iii) identify the vital environmental parameters which influence habitat prediction and, (iv) recommend conservation strategies based on the area suitability of these species.

Materials and Methods Materials and Methods Materials and Methods Materials and Methods
Description of the study site Study area northern India, Uttar Pradesh (UP) (Figure 2) is situated between 23°52'N and 30°24'N latitude and 77°5' and 84°38'E longitude. Its geographical area of 240928 km 2 is divided into four ecozones: Tarai, Bundelkhand-Vindhyan, Semi-arid and Gangetic with varying combination of vegetation, temperature and precipitation (Islam and Rahmani, 2004;Jha, 2015). It has a humid subtropical dry winter climate (FSI, 2017). The ranges of elevation, mean annual temperature and mean annual precipitation are 23-935 m, 21.01 o C -26.56 o C and 601-2145 mm, respectively (Fick and Hijmans, 2017). The winter conditions here prove to be favourable for the wintering vultures which look to escape the harsh winter conditions in their breeding grounds. UP has tropical moist and dry deciduous forests as well as tropical thorn forests, hilly terrain with cliffs, agriculture landscape and abundant water bodies, all considered as important components of suitable vulture habitat. In the past two decades sightings of EBV and EGV have been reported from the state during the winter months of November to March (Jha, 2015;Sullivan et al., 2009;iNaturalist users and Ueda, 2020). in different eco-zones. Data sources: Fick et al. (2017) and Buchhorn et al. (2020) Model selection Stockwell and Peterson (2002) suggested that 50-100 data points should be reasonable for robust model prediction. In practicality many times this sample size is not available especially in the case of rare and threatened species, hence, Aguillar et al. (2016) suggested implementation of various algorithms. In the present study EBV as well as EGV both had much lower number of occurrence points, therefore, all the three categories of SDM modelling tools were selected, for example, regression models (GBM, GLM, MARS), classification techniques (CTA), and machine learning (ANN, MaxEnt, RF, SVM) as recognised by some researchers (Bucklin et al., 2015;Zang and Li, 2017;Fruh et al., 2018;Passadore et al., 2018 etc.). All these SDM models work on the principle of likelihood of animal occurrences at a spatial location using environmental variables, thus quantifying the environmental conditions that may lead to species occurrence (Hirzel and Le Lay, 2008;Farrell et al., 2019). Therefore, the models required vulture occurrence points (longitude and latitude) and environmental or bioclimatic characteristics of each occurrence points as inputs.

Model inputs
Occurrence (presence) data for both EBV and EGV was collected and compiled from four sources namely, primary data collected from (i) field surveys along the transects (Figure 3) in all ecozones covering all important landcovers, and secondary data from (i) published literature (Jha, 2015), (ii) official government records {census records from selected forest divisions and whole UP state (UPFD and BNHS, 2021)}, and (iii) research grade citizen science records from eBird (www.ebird.org) (Sullivan et al., 2009) and iNaturalist (www.inaturalist.org; iNaturalist users and Ueda, 2020). As the data was compiled from different sources, expert as well as non-expert, it was prone to some biases such as recurrent sightings, crowded observations and latitudinal bias. These were corrected with the help of the duplicate removal tool of MS Excel, Spatial rarefication tool of the SDM tool box (Brown et al., 2017) and creating a bias file. The rarefying distance was set at 4km according to the shelter sites defined by . The final rarefied data set had 16 points for EBV and 17 points for EGV from 20 unique points each. Dashed lines indicate railway transects The second set of input was the environmental variables as raster layers {19 bioclimatic variables (Table  1) downloaded from www.worldclim.org (Fick and Hijmans, 2017), elevation data Earth Resources Observation and Science (EROS) Center, 2017, two sets of Normalized Differentiated Vegetation Index (NDVI) data of winter (January) and summer (June) seasons (Didan, 2015) download from www.earthexplorer.usgs.gov, and land use and land cover (LULC) downloaded from www.land.copernicus.eu (Buchhorn et al., 2020)}. LULC types were merged and the map was reclassified into six classes (Forest, Water, Scrubland, Agriculture, Built-Up and Wasteland) as the present study did not require a detailed classification (Jha and Jha, 2021). . .
. The four environmental layers were chosen to cover the shelter and foraging components of the niche. All 23 variables were projected to the WGS 1984 Geographical Coordinate System (GCS) and then clipped using a shapefile of the political borders of UP. Finally, all the layers were rescaled to 1 km (30 sec) spatial resolution. To remove the collinearity among the variables Pearson's correlation test was carried out using the

Habitat classification and mapping
Each model generated a heatmap which showed a continuous scale of habitat suitability from 0 (unsuitable) to 1 (highly suitable). This continuous scale was reclassified into four categories: unsuitable (0.0-0.25), low suitability (0.25-0.50), moderate suitability (0.50-0.75) and high suitability (0.75-1.00) as suggested for raptors and vultures (Zhang et al., 2019;Jha and Jha, 2021a). With the help of this categorisation habitat suitability maps were prepared using ArcGIS 10.3 and suitability area was calculated using pixel number and value.

Results and Discussion
Results and Discussion Results and Discussion Results and Discussion

Model performance
In total eighteen models were developed, nine each for EBV and EGV. The results of model evaluation in terms of AUC, Kappa and TSS for each species are presented in Table 2. A corresponding performance matrix is also presented in Figure 5. The value of model evaluators for a species varied a lot between models.
Such changes in performance have been reported earlier also (Grenouillet et al., 2011). The AUC values for EBV ranged from 0.7 to 1.0 among nine model tools tried while for EGV it ranged from 0.6-0.8. As per the AUC values, the models fell in the excellent, good and fair categories for EBV and good, fair and poor for EGV. AUC is considered to be a reliable indicator of model robustness and has been used in many studies for the same (Malabet et al., 2020;Abdelaal et al., 2019;Abolmaali et al., 2018). However, some studies have also indicated that AUC values may be misleading and may denigrate the accuracy of the model (Fourcade et al., 2014;Jiménez-Valverde, 2014). Therefore, alongside AUC, the additional model evaluators like, Kappa and TSS values were also recommended and used to judge the performance of the models (Allouche et al., 2006;Shabani et al., 2018;Zhang et al., 2020). In the present study, the Kappa values for the EBV ranged from 0.0 to 1.0 but for EGV, this ranged from 0.0 to 0.4. The TSS value for EBV ranged from 0.4 to 1.0 while for EGV the range was from 0.2 to 0.6. As per the Kappa values the models were in the excellent, good and poor categories for EBV and fair and poor for EGV. As per the TSS values the models fell in the categories of excellent, good and fair for EBV and good, fair and poor for EGV. In general, EBV showed better degree of model performances or robustness than EGV. However, some of the models fell in good enough categories (fair and above) to be relied upon in both the species (Figure 5). Table  Table Table  Table 2

Model selection
High performance of model evaluators (Fair, good and excellent) produced the candidate models for further selection for both the species. Such selection was based on minimum two performers out of three (AUC, Kappa, TSS). Thus, such candidate models for EBV were ANN, CTA, GBM, GLM, MARS, MaxEnt, RF and SVM while for EGV they were GLM, GBM and MARS. This indicated that for the former species all the three types of algorithms (regression, classification and machine) were good but for the latter species only regression algorithm proved to be good. This is contrary to the general notion that machine learning algorithms/complex models, being advanced, should perform better than regression algorithms (Li and Wang, 2009). However, further refining of this selection was needed since a lot of variation in prediction of suitable area among the models was observed. Ground-truthing of the output, apart from model evaluators, formed another basis of selection or rejection of models. In this regard, the heatmaps (Figure 6 and 7) generated by MARS for both the species were found to acutely under predict the existing distribution of both EBV and EGV. Therefore, MARS was also rejected, even though it was found useful as per model evaluators. Such observations and corrections have been recommended by Mi et al. (2017) who noted that modelers should not depend on some model evaluating metric alone, but also should base their assessments on the combined use of visualization and expert knowledge. Thus, it is hypothesised that barring MARS all other models (ANN, CTA, GBM, GLM, MaxEnt, RF and SVM) were found to be acceptable for EBV and only GBM and GLM were suitable for EGV in the conditions of the present study.

Ensemble and component models
The ensemble models of EBV ('Good' to 'Excellent') and EG (all 'Fair') both had relatively better range of values in all the three categories of model evaluators in comparison to individual models ( Figure 5). The use of ensemble modelling methods for species distribution modelling have been promoted over single algorithms since, using the ensemble approach is thought to reduce model uncertainty and increase its robustness in modelling species distributions accurately (Thuiller, 2003;Araújo and New, 2007;Marmion et al., 2009;Kaky et al., 2020). In ensemble approach the results of individual algorithms are combined (Passadore et al., 2018), however, the combination of all available models is not advised, because there will always be some with weaker or instable output, which disqualify the output of the whole ensemble (Fruh et al., 2018). Present findings of ensemble include only those algorithms which performed Fair, Good and Excellent as per AUC. Thus, the ensemble of EBV contained ANN, CTA, GBM, GLM, MARS, MaxEnt, RF and SVM algorithms while of EGV contained ANN, CTA, GBM, GLM, MARS and RF. Former species included all the models but in the case of latter two machine learning algorithms, MaxEnt and SVM were excluded. Generally, machine learning tools perform better than simple regression tools but contrasting results are also reported since model performance is known to depend on many factors like, sample size, locality factors, spatial extent of the study area, number of ecologically and statistically significant variables, interaction of species with environment etc. .

Habitat suitability
The total area available for both EBV and EGV in UP is 240928 km 2 but large parts of this were found to be unsuitable habitat. The available area according to different categories of suitability varied as per predictions by different models for both the species. The detailed and comparative expanse of these areas is presented species-and model-wise in Table 3 and Figures 6 and 7. The suitable area for EBV ranged from a minimum of 1% (MARS) to a maximum of 36% (SVM). The average unsuitable area and suitable area of all prediction models were found to be 82% and 18%, respectively ( Figure 6 and Table 3). The suitable area for EGV ranged from a minimum of 1% (MARS) to a maximum of 67% (RF). The average unsuitable area and suitable area of all prediction models were found to be 72% and 28%, respectively (Figure 7 and Table 3). A key observation of this study was that the expanse of area, as predicted by the different models, as well as the values of different model evaluators varied considerably from one model to another. Due to this no consensus could be reached as to which model was the best among all acceptable models. Therefore, the ensemble of all superior models was chosen as the most suitable since it combined the positives of all the component models (Stohlgren et al., 2010;Thuiller and Münkemüller, 2010;Aguilar et al., 2016;Fruh et al., 2018).
For these two species (EBV and EGV), the ensemble was made with all the superior component models, mentioned in previous section except MaxEnt in EBV. MaxEnt was excluded as it was found that at default settings, the model gave almost equal weightage to all the habitat determining variables which was not realistic. Models that failed the test of ground truthing (i.e., MARS) were also included in the ensemble as their model evaluators indicated it a strong model. Therefore, it was surmised that MARS was a conservative model for habitat prediction.
The ensemble model predicted 79.6% area as unsuitable (191758 km 2 ), 9.9% as area with low suitability (23845 km 2 ), 8.3% area with moderate suitability (20171 km 2 ) and 2.1% area with high suitability (5153 km 2 ) for EBV (Figure 8; Table 3). Similarly, for EGV, the ensemble model predicted 84.5% area as unsuitable (203670 km 2 ), 14.1% as area with low suitability (33922 km 2 ), 1.4% area with moderate suitability (3334 km 2 ) and very small area with high suitability (2 km 2 ). Among the two species, EBV had more suitable habitat for wintering as compared to EGV. However, since UP is a large state, 20.3 % (49169 km 2 for EBV) and 15.5% (37258 km 2 for EBV) suitable area is sizable for conservation if coupled with good management practices.
The ensemble habitat suitability map when overlaid with a map of the ecozones of UP (Jha, 2015) also provided useful insights about the habitat of both species. The Tarai ecozone proved to be the most favourable for both the species, followed by the Vindhyan-Bundelkhand, Gangetic and Semi-Arid ecozones, respectively. The Tarai ecozone is characterised by tropical moist deciduous forest and a combination of lower temperature and higher precipitation as compared to the rest of the state. This indicates a suitability of this combination for the habitat of these two species. When compared to the land use land cover map, it was also evident that these species overwhelmingly preferred forests for their habitat over agricultural areas, with very little suitable area being found in the agrarian part of ecozones, primarily the Gangetic ecozone. This also suggested shyness of these vultures around human presence which is in stark contrast to another vulture species found in UP, the Egyptian Vulture, as the latter has been recorded on riverbanks, garbage dumps and electricity poles (Bahadur, 2019).
However, it is proposed that the highly suitable area where speculated good conditions for shelter, foraging and water requirement are available, must be given high priority of protection. In general, a decline of Eurasian vultures has been blamed on increased agricultural activities, the decline of wild ungulates, persecution etc. (Campbell, 2015). Therefore, no habitat degradation and reduction in ungulate population should be allowed in highly suitable areas. After thorough investigation, moderate and low suitability areas, wherever necessary, may be enriched with shelter and forage needs, by raising tall tree plantation and rearing livestock. The unsuitable areas may also be augmented with food supply to provide a safety-net to these vultures. Table  Table Table  Table 3 3 3 3. . . . Species-wise habitat suitability area (km 2 ) in four classes predicted by different SDM algorithms at default settings  The Pearson test resulted in elimination of 11 variables (bio4, bio5, bio6, bio7, bio10, bio11, bio12, bio16, bio17, bio18, and bio19) from model prediction due to collinearity. The percentage contribution to model building for all non-corelated 12 variables was available as an output of each model. The five most important variables (highest cumulative contributors) for each species along with their percentage contribution are given in Figures 9 and 10. These important variables contributed more than 50% (52-95%) in model prediction except CTA in EBV and MaxEnt in both EBV and EGV. Variable prediction in the case of CTA for EBV (only two variables) and ANN for EGV (only three variables) appeared to be an outlier since these few variables contributed 100%. But this aberration was not recorded in the Ensemble models of either species.
Identification of important variables as vital contributors have been emphasized in earlier research on vultures in Central India and Gangetic-Deccan-Thar region of India Jha, 2020, 2021a,b) and elsewhere (Phipps et al., 2017;Zhang et al., 2019;Anoop et al., 2020). In the ensemble model which is the consensus of superior models, five top variables played their role in model determination up to 78% (EBV) and 57% (EGV) (Figure 11). Among these five, NDVI was the most important environmental variable for both EBV and EGV. Since NDVI can be used as a proxy for food source for vultures (Santangeli et al., 2018), it is one of the most important concerns for any non-breeding species. The most important variable of precipitation for both species was bio13 (precipitation of the wettest month). While the wettest month in UP (July/August) does not coincide with the visit of these species, it has an important role to play in determining the habitat in the winter season. The rainfall received during the monsoon period determines the moisture content of the vegetation and while the rivers are perennial and snow-fed, the water content, for the rest of the year, of the other water bodies such as lakes and reservoirs depends on the monsoon. Water is also an important habitat determinant for vulture species as vultures prefer feeding and roosting sites in close proximity to water-bodies Kanaujia and Kushwaha, 2014). Finally, bio09 (mean temperature of the driest quarter) and bio08 (mean temperature of the wettest quarter) were the most important variables of temperature for EBV and EGV, respectively. However, the influence of bioclimatic variables in habitat determination, in general, has been noted in other studies on raptors and scavengers (Gschweng et al., 2012;Liminana et al., 2012;Phipps et al., 2017;Zhang et al., 2019;Anoop et al., 2020;Jha and Jha, 2021a).     Interestingly, the results on the important role of variables in the present study does not support an earlier hypothesis that LULC is the most important contributor in model prediction in central India in the case of migratory vultures (EBV, EGV and Himalayan Griffon) as derived by MaxEnt algorithm . Jha and Jha (2021a) also had similar findings as of , contrasting to the present study, when all the three migratory vultures were combined together. It is speculated that the modelling tool selected and its settings could be the main reason and the ensemble is not an exact reflectance of MaxEnt or any of its component models. However, to ascertain this, further research with only MaxEnt could be attempted.

Conclusion Conclusion Conclusion Conclusion
This is among very few habitats suitability studies in India and the first one using ensemble modelling and comparing the performance of different models for Eurasian vultures. This study provides information on suitable models and habitat details of wintering Eurasian vultures in northern India, UP which can aid in vulture management and conservation. The superior species distribution models for EBV were SVM, ANN, GBM and RF while GLM and GBM for EGV. However, ensemble modelling provided the better option for habitat projection for both EBV and EGV in UP. Tarai ecozone was most suited for both the vultures followed by Vindhyan-Bundelkhand, Gangetic and Semi-Arid ecozones. The highly suitable areas concentrated in the northern part of the state which can be designated as vulture protection area and be prioritised for stopping any kind of degradation and reduction in herbivore population. Medium and low suitability areas may be enriched with shelter and forage needs, by raising tree plantation and rearing livestock. The unsuitable areas may also be augmented with food supply to provide safety-net to these vultures.