Modeling global habitat suitability and environmental predictor of distribution of a Near Threatened avian scavenger at a high spatial resolution

Vultures are among the most vulnerable birds in the world. The bearded vulture (Gypaetus barbatus) is among the threatened species of vultures and listed as Near Threatened. The species is widely distributed across the Palearctic, Afrotropical, and Indomalayan regions. The species faces several threats such as poisoning, direct persecution, habitat degradation, and collisions with powerlines and wind power farms. Thus, knowing the global habitat suitability of the species and environmental predictors of the species distribution can facilitate the species conservation. In this study, we applied a maximum entropy approach, 10,585 distribution records, and 10 environmental variables to model the bearded vulture's global habitat suitability at high spatial resolution [30-arc-second (1 km)]. We also estimated protected area coverage for the species' suitable habitats. We identified 8,117,231 km2 of suitable habitat for the species across its global range in Europe, Asia, and Africa. The results showed that topographic diversity is the most important predictor of the species distribution across its distribution range. Results of estimating the area of suitable habitats of the bearded vulture within protected areas revealed that only 16.26% of the species' suitable habitats are protected. The areas that were identified to have the highest suitability for the species have high priority for the conservation of this iconic species thus these areas should be included in the network of protected areas.


. Introduction
Human activities over the past decades have led to habitat destruction and a decrease in numerous species (Hanski, 2011;Hansen et al., 2013;Haddad et al., 2015). Climate change, land use change, and road development are among the human activities that are causing the decline in biodiversity worldwide (Thomas et al., 2004;Hansen et al., 2013;Haddad et al., 2015;Newbold, 2018). Among the different groups of species, some are more sensitive to human activities such as raptors (hawks, harriers, kites, eagles, falcons, owls, and vultures). This group provides critical ecosystem services such as disease mitigation, agricultural production, and waste-disposal services (O'Bryan et al., 2018;Carucci et al., 2022). Among raptors, vultures are the most threatened and sensitive group (Botha et al., 2017). In fact, due to human persecution, inadvertent poisoning, and mortality associated with human infrastructures [i.e., roads, buildings, and wind power farms (Ogada et al., 2012;Safford et al., 2019)], vultures have become one of the most threatened taxa on the earth (Botha et al., 2017).
The bearded vulture (Gypaetus barbatus), a cliff-nesting bird of prey, is among the most threatened species of vultures and is listed as Near Threatened by the IUCN (BirdLife International, 2022). The species is widely distributed across the Palearctic, Afrotropical, and Indomalayan regions, with an area of ∼61,700,000 km 2 (BirdLife International, 2022). There has been a major decline in the distribution range of the habitat specialists' avian scavengers such as bearded vulture, Gypaetus barbatus (Acharya et al., 2010;Safford et al., 2019). As one of the most well-known vulture species, the bearded vulture has experienced a severe population decline in many of its distribution regions and may soon become endangered or reach the status of local extinction Sheykhi Ilanloo et al., 2020). However, population trends varied significantly throughout the species range; while the European population species has increased since 1980, some other local populations experienced a decline of 89.3% in Nepal (BirdLife International, 2022). At various phases of their life, bearded vulture chooses altitude grasslands and steep pastures as their primary habitat. When looking for food, bearded vultures often soar near mountain slopes or over the ground in their visual search for food, looking for bones and other remnants of carcasses (Margalida, 2008). Only free landscapes allow for this type of search for food (BirdLife International, 2022). Since the bearded vulture represents a significant umbrella species, any measures taken to ensure the survival of this species will also benefit other species that co-exist with it. As one of the most charismatic vulture types (Aguilera-Alcalí et al., 2020), the bearded vulture can only be found in the highest mountain ranges in Asia, Europe, and Africa (Orta et al., 2020). Our knowledge of the bearded vulture distribution varies from one mountain range to another one across the species distribution range (BirdLife International, 2022). While the species distribution is well documented in the Central European Highlands and Southeast European Highlands, it is not well known in other mountain ranges across its global distribution ranges, particularly in Africa (BirdLife International, 2022).
This study aimed to model the global geographic distribution of the bearded vulture, to determine the most important environmental factors affecting the global habitat suitability, and to estimate protected area coverage for the suitable habitat of the species. It is known that at large scales, climate serves as the most important factor in shaping species distribution (Pearson and Dawson, 2003). Thus, we are expecting climatic factors to be the most influential predictor of the bearded vulture's global distribution. Human activities during recent decades are a further factor affecting the species distribution ranges (Newbold, 2018;Xu et al., 2019). For instance, environmental changes triggered by humans via their land usage and the resulting habitat degradation caused a severe decline in the distribution ranges of several species. Hence, we are expecting the human footprint index to play a key role in species distribution. The results of this study will facilitate large-scale conservation planning of this iconic and rapidly declining avian scavenger.
. Materials and methods

. . Species global presence data
Global distribution records (Figure 1) of the bearded vulture ( Figure 2) were collected from online databases such as the Global Biodiversity Information Facility (GBIF, 2022), eBird (https://ebird.org/home), and VertNet (http://vertnet.org/). In total, 43,670 distribution records were collected, and then duplicates were deleted. After removing duplicates, 17,506 presence points remained. We then assessed the representation of the species distribution records within major mountain ranges across the species' global distribution range (Snethlage et al., 2022).

. . Environmental variables
To map the global habitat suitability of the bearded vulture, we considered 14 environmental variables (Table 1) related to climate, topography, vegetation, and human presence. Climatic variables were obtained from the Wordclim at 1 km spatial resolution (Fick and Hijmans, 2017). To consider topography in our model, we estimated topographic heterogeneity based on the Shuttle Radar Topography Mission (SRTM) elevation model (Jarvis et al., 2008). For vegetation, we used the normalized difference vegetation index (NDVI). The human footprint index was used to quantify the anthropogenic impact on the habitat suitability of the species (Venter et al., 2016a,b). This index was created by combining data on the built environments, population density, electric infrastructure, crop lands, pasture lands, roads, railways, and navigable waterways (Venter et al., 2016a,b). We performed a VIF test and only included variables (10 variables) .
/fevo. .  with VIF values less than 10 in the Maxent model (Guisan et al., 2017). Considering the spatial resolution of environmental variables (1 km resolution), we thinned the dataset to ensure that observations were at least 1 km apart. Thus, 10,585 distribution records with at least a 1 km distance were used in habitat suitability modeling.

. . Global habitat suitability modeling
Maximum entropy modeling (Maxent) is the most popular algorithm among many algorithms for modeling species habitat suitability. Like the other habitat suitability algorithms, Maxent uses presence data from the target species and predictors that .

Variables Abbreviations
Mean diurnal range Bio2 Isothermality Bio3 Temperature seasonality Bio4 The mean temperature of the wettest quarter Bio8 The mean temperature of the warmest quarter Bio10 Precipitation of the driest month Bio14 Precipitation of the coldest quarter Bio19 Topographic heterogeneity TH

Human footprint HF
Normalized difference vegetation index NDVI present environmental conditions to find suitable habitats for the target species across the study area. In this study, the Kuenm R package (Cobos et al., 2019) was used for modeling the bearded vulture's global habitat suitability. We applied this package to create Maxent candidate models with multiple combinations of regularization multipliers, feature classes, and sets of variables. Then, the best parameters for modeling were selected based on the statistical significance, predictive power, and model complexity (Cobos et al., 2019). The area under the ROC curve (AUC) was used to evaluate the species distribution model's performance (Swets, 1988

. . Protected area coverage for the bearded vulture
To estimate the protected area coverage of suitable habitats of the bearded vulture across its global range, the species continuous habitat suitability map was converted into a binary suitableunsuitable map (Guisan et al., 2017). The 10-percentile training presence threshold was used to convert the continuous map into binary (Phillips et al., 2006). Then, the binary habitat suitability model was overlaid on the protected area layer. Protected areas' data were obtained from the protected planet (www.protectedplanet.net) (UNEP-WCMC and IUCN, 2022).

. . Distribution in the major mountain range
We showed that the breaded vulture occurs in 26 major mountain ranges. Table 2 shows the number of distribution records in each major mountain range across the species distribution range. The highest number of distribution points was recorded in the Central European Highlands.

. . Global habitat suitability model
According to the AUS value (AUC = 0.924), the global habitat suitability model developed for the breaded vulture (Gypaetus barbatus) shows high performance. The Maxent model shows areas with high suitability for the species across Europe, Asia, and Africa. The model is created based on distribution records of the species collected from 2000 to 2022 across the species distribution range. We estimated 8,117,231 km 2 of suitable habitat for the species across its global range in Europe, Asia, and Africa (Figure 3).

. . Variable importance
Results of estimating the importance of the environmental variables to the bearded vulture global habitat suitability model . /fevo. . showed that topographic heterogeneity was the most important variable (with 46.2% contribution) in predicting suitable habitats of the species with a positive relationship (Table 3). Areas with high topographic heterogeneity have higher suitability for the species. Precipitation of the driest month (Bio14) was the second most important variable (with 21.6% contribution) in shaping the species' habitat suitability with a positive influence meaning that areas that receive more precipitation during the driest month have high suitability for the species. Figure 4 shows the results of the jackknife test of variable importance. The environmental variable with the highest gain, when used in isolation, is topographic heterogeneity, which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is Bio14, which therefore appears to have the most useful information that is not present in the other variables (Figure 3).

. . Representation of suitable habitats within protected areas
Results of estimating the area of suitable habitats of the bearded vulture within protected areas showed that 1,320,533 km 2 of suitable habitats of the species are protected. Considering that area of suitable habitat for the species is 8,117,231 km 2 , thus 16.26% of the suitable habitat of the bearded vulture is currently protected ( Figure 5).

. Discussion
Mountain ecosystems are hotspots of biodiversity and have high priority for conservation (Antonelli et al., 2018;Rahbek et al., 2019). These ecosystems are home to many iconic species in need of urgent conservation actions like the bearded vulture. This study presents a high-resolution global scale habitat suitability map for the bearded vulture to help large-scale conservation planning of the species. The overall extent of the occurrence of the species . /fevo. .

FIGURE
Results of the jackknife test for variable importance. is estimated to be 61,700,000 km 2 (BirdLife International, 2022). However, our model on the global distribution of the species showed that the area of suitable habitat for the species is much lower than what was previously thought (BirdLife International, 2022). This highlights the importance and usefulness of SDMs for precise estimation of the area of suitable habitat for species with conservation importance. Despite we were expecting that climate will outperform other variables in predicting the species' global habitat suitability, the results showed that topography was the most influential factor. Topographic heterogeneity is an indicator of habitat diversity and species diversity (Stein et al., 2014;Badgley et al., 2017), which, in turn, influences food resources availability to bearded vultures (Sheykhi Ilanloo et al., 2020). Vignali et al. (2021) built several habitat suitability maps with respect to bearded vulture age class (immature and adult) and season (warm and cold) in the Swiss Alpine range. In most of the models, the environmental predictor that contributed most to explaining the habitat suitability of the bearded vulture was food availability. The precipitation of the driest month was the second most important predictor of the species distribution range showing the importance of precipitation for the species' habitat suitability. Sheykhi Ilanloo et al. (2020) .
/fevo. . modeled the habitat suitability of the bearded vulture in the Kopet Dagh Mountains in Asia and found that the precipitation was the most influential predictor of the species' suitable habitat in the area. Surprisingly, results of variable importance revealed that the human footprint index did not play an important role in shaping the species distribution showing that at least at a global scale, the species' suitable habitat is not affected by human activities. Confirming previous findings, our results show that species responses to environmental factors such as climate, topography, vegetation, and anthropogenic factors are speciesspecific and depend on the scale of studies (Roland and Taylor, 1997;Holland et al., 2004;McGarigal et al., 2016). We showed that while topography shaped the species' global distribution range, precipitation was the most important determinant of the species distribution in the Kopet Dagh Mountains in Asia (Sheykhi Ilanloo et al., 2020), and food availability was the most important variable in explaining habitat suitability of bearded vulture in the Swiss Alpine range. In comparison with other vulture species, Panthi et al. (2021) modeled the global habitat suitability of the Egyptian vulture (Neophron percnopterus) and found that climate is the most important determinant of the species' global habitat suitability. Protected areas are important tools for the conservation of biodiversity (Gaston et al., 2008;Mawdsley et al., 2009;Watson et al., 2014;Gray et al., 2016). Threats that species generally face outside protected areas are absent or very low in protected areas (Leberger et al., 2020). Thus, knowing to what extent species' suitable habitats are located within protected areas can help species conservation and future selection of protected areas. In this regard, SDMs are very useful tools to investigate the effectiveness of protected areas in the conservation of species and identify gaps in protected area coverage for species' suitable habitats (Mudereri et al., 2021;Escobar-Luján et al., 2022). Here, we applied SDMs to estimate protected area coverage for the suitable habitat of the species and found that <17% of the bearded vulture's suitable habitat is currently protected. This highlighted the necessity of new protected area establishment in mountain ecosystems in general and suitable habitats of this iconic and threatened avian species in particular.
Vultures are among the most threatened birds that are in need of urgent conservation actions (Botha et al., 2017). Like many other vultures, the bearded vulture has been listed as Near Threatened because it has undergone a moderately rapid population decline over the past three generations (BirdLife International, 2022). We found a high number of distribution records for the species in the following mountain ranges: Central European Highlands, Southeast European Highlands, Himalayas, Caucasus Mountains, and Anatolian/Armenian Highlands showing that more efforts have been granted for monitoring of this species or avian species in general in these mountains. While few observations have been reported in other mountain ranges such as Iranian Plateau, that the species is resident and common there (Sheykhi Ilanloo et al., 2020). Thus, more effort should be granted in data poor mountain ranges to monitor the species populations, particularly by using GPS transmitters and record their presence, which is a critical step in understanding the species' global distribution and conservation (Di Vittorio et al., 2018). Moreover, as suggested by Di Vittorio et al. (2018), the development of citizen science programs can significantly increase our knowledge about the species distribution and population in data poor regions. We particularly encourage awareness-raising programs to increase local people's awareness about the species and the species' highly suitable habitats identified in this study. The species population in Europe has increased due to conservation actions, thus aforementioned programs and conservation actions are necessary for the highly suitable habitat of the species in Asia and Africa in which the species experiencing dramatic population decline (BirdLife International, 2022). Species distribution models have been frequently used in studying avian ecology and conservation (Engler et al., 2017;Fourcade et al., 2017;Brambilla et al., 2018Brambilla et al., , 2019Brambilla et al., , 2020Ramellini et al., 2019;Burns et al., 2020;Ferrer-Paris and Sánchez-Mercado, 2021;Li et al., 2021;Lu et al., 2021;Mudereri et al., 2021;Avotins et al., 2022;Condro et al., 2022;Escobar-Luján et al., 2022). However, studies that applied SDMs in studying avian global distribution are still limited (Engler et al., 2017;Ferrer-Paris and Sánchez-Mercado, 2021;Stiels et al., 2021). In this study, we modeled the global habitat suitability of the bearded vulture and showed that the species distribution range is very smaller than what was previously thought. In addition, we showed that a small proportion of the species' suitable habitats is located within protected areas. We believe that the establishment of new wind power farms should be avoided in highly suitable habitats or should be established with caution to reduce the species' collision with turbines (Carrete et al., 2012;Murgatroyd et al., 2021). Areas that are identified as the most suitable habitats for the species should be prioritized for the conservation of this charismatic vulture species.

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions
MY conceived the ideas, designed the methodology, collected the data, and wrote the original draft with inputs from SM (Section 1). AK conducted the analyses. MY, SM, and AK reviewed, edited the text, and gave their final approval for submission. All authors contributed to the article and approved the submitted version.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.