Modelling the distribution of bat activity areas for conservation in a Mediterranean mountain range

There is a lack of studies designed to detect the most important areas for bat conservation. In this context, areas of high bat activity have been rarely considered in the delimitation of protected areas for bats, which are generally focused on the protection of roosting sites. This has been due to the difficulties of sampling the distribution of these nocturnal animals when moving at night. This methodological constraint has been overcome by the development of bioacoustic sampling, which allows mapping the occurrence of active bats over large areas. In this study, we use bat detectors to sample the distribution of bat activity in central Spain. This region is under the environmental effects of a mountain range (Guadarrama Mountains) and the urban encroachment of the city of Madrid. The occurrences provided by the detectors were used to produce species distribution models of which the resulting layers were arranged to detect the most suitable areas for bat richness and rarity indices. We performed a gap analysis to explore whether the areas most commonly used by active bats are covered by the current network of protected areas. The results showed that the best areas of high bat activity are located at the piedmont of the mountains and that most of these areas overlap with the existing network of protected areas. The best areas for bats excluded the most urbanized areas and within a similar urban gradient, protected areas tended to be located within the best sites for conservation. These results suggest that bats currently benefit from a network of protected areas initially aimed to protect birds and habitats (Natura 2000). In addition, monitoring areas of high bat activity could complement roosting site protection in the conservation of bat assemblages.


Introduction
Bats comprise the second most species rich order of mammals with approximately 1,400 species distributed throughout the planet (Zukal, 2020). This taxonomic diversity makes the group an interesting target of conservation, especially because they are increasingly affected by the pervasive effect of habitat loss, pesticides, wind farms and the alteration of roosting areas (Mickleburgh et al., 2002;Voigt and Kingston, 2016). In this context, the location and protection of suitable areas for bats (e.g., both roosting sites and foraging areas) have been considered a suitable approach to their conservation (Russo and Jones, 2003;Medell ın et al., 2018). The location and protection of bat roosts is a common conservation policy in many countries (EUROBATS, 2006) but, unlike other, easier to detect taxonomic groups (e.g., birds, butterflies and plants), there is a lack of cartographic approaches devoted to detecting the most important habitats and/or areas for the conservation of bats (Razgour et al., 2016). Since the foraging areas can be located far away from resting sites, it can be assumed that bats will track food resources over large areas (Medell ın et al., 2018;Nad'o et al., 2019) and that this ubiquity will present an opportunity to map the most suitable areas for bat activity.
Locating areas of high bat activity has been difficult due to challenges of detecting active bats at night. However, this restriction has been overcome by the improvement of bioacoustic approaches for bat detection (Walters et al., 2012). Bat detectors allow sampling of active bats over large areas and the use of resulting occurrences to produce species distribution models (Elith and Leathwick, 2009;Razgour et al., 2016). In this way, the resulting maps can be used to explore the distribution of the most suitable areas for bat activity.
Here, we study bat distribution in central Spain, an area located within the highly biodiverse Mediterranean region (Myers et al., 2000) where bat distribution is still poorly known (Palomo et al., 2007;De Paz et al., 2015). The area is crossed by a mountain range (Guadarrama Mountains), which induces an upwards decrease in temperature and increase in precipitation from the surrounding lowlands (Gonzalez-Hidalgo et al., 2016). This climate turnover produces an altitudinal succession of vegetation belts and a concomitant change of animal assemblages (e.g., Ruiz-Labourdette et al., 2012;Flores et al., 2018). In a context of global change in which the Mediterranean is under the pervasive effects of increasing temperature and decreasing precipitation (Giorgi and Lionello, 2008), the region is an interesting setting in which to assess the potential role of protected areas to conserve bats in the future. In addition, the region is under the effect of the metropolitan area of Madrid, experiencing considerable urban encroachment (Hewitt and Escobar, 2011;Fig. 1). It can be thus conjectured that bat distribution will also be stressed by this process of urban encroachment (Goddard et al., 2010). Therefore, we study bat distribution using four complementary approaches: a) we sample the distribution of active bats within the study area to produce species distribution models (Elith and Leathwick, 2009); b) we use the occurrence probabilities of the species to detect the most important areas for bat conservation (Razgour et al., 2016); c) we carry out a gap analysis (Scott et al., 1993) to detect whether the best areas for bat conservation are included in the regional network of protected areas (Buckman-Sewald et al., 2014;Bosso et al., 2016;Kerbiriou et al., 2018); and d) we explore whether bats avoid the most urbanized areas (Tena et al., 2020a) and whether protected areas prevent the potential effects of urban encroachment.

Study area
The study area is located in central Spain, in a region of around 10,000 km 2 divided from NE to SW by the Guadarrama Mountains (Fig. 1). These mountains, which range from 600 to 2400 m a.s.l., are covered by an altitudinal succession of cereal fields, gum rock shrublands Cistus ladanifer and schlerophylous Holm oak woodlands Quercus ilex in the piedmont (under 1000 m) to Scots pine forests Pinus sylvestris and broom shrublands Cytisus oromediterraneus in the upper parts of the mountains (above 1400 m). Between 1000 and 1400 m, the vegetation is composed of shrublands dominated by laurel-leaf cistus Cistus laurifolius interspersed with Pyrenean oaks Quercus pyrenaica. Common ash Fraxinus excelsior woodlots interspersed with meadows occur in wet valleys within this mountain range (Teller ıa, 2019). Most of the habitats other than cereal fields are managed as pasturelands for extensive cattle breeding except in the case of pinewoods, which are mainly managed for timber production. Most highlands (over 1700 m a.s.l) of the Guadarrama Mountains were designated a National Park in 2013 (L opez and Pardo, 2018). The surrounding piedmont has been covered by an extensive network of buffer areas and reserves resulting from different EU directives (e.g., Natura 2000, Fig. 1). A large part of the region is densely populated by humans (6.5 million inhabitants; Cincotta et al., 2000) and occupied by a dense network of infrastructure (roads, railways, power lines) and residential areas. As a result, a marked NW-SE gradient is defined between the less urbanized mountain areas and lowland areas with high human population density around the city of Madrid (Fig. 1).

Bat sampling
Bats were sampled across 274 sampling points distributed within the study elevation gradient (Fig. 1) during the breeding season (June-August) of the years 2014 (58 sampling points), 2015 (82) and 2016 (134). In this sampling we tried to track the environmental variability of the studied area, which led us to select a sample of more elevated and treecovered location than the average according to a regional random sampling (Appendix S1). Bats were recorded by different ultrasound bat detectors (Echo Meter 3, Song Meter 2 Wildlife Acoustics). Each location was visited four times (twice in July, once in June and once in August) within the first three hours after sunset (a period of high bat activity; Barataud, 2012). Finding just one sequence during any of the four visits was enough to assume the presence of a species at that sampling point. All sequences were recorded as full-spectrum in WAV format and filtered to detect bat calls using the program Kaleidoscope (Wildlife Acoustics, Inc.). The filter settings were specified to detect any signal between 8 and 120 kHz, from 2 to 500 ms and with a minimum of two calls per sequence. Batch function split each sequence in a maximum duration of 5 seconds. We then analysed the WAV files using BatSound 4 (Pettersson Elektronik AB, Uppsala). The recordings were analysed using a sampling frequency of 44.1 kHz, with 16 bits/ sample and a 512 pt. Fast Fourier Transform with a Hamming window for analysis. At least two random echolocation calls were analysed manually from each sequence to identify the species (Rydell et al., 2017). To do this, we measured the following parameters (Russo and Jones, 2002;Barataud, 2012): start frequency, end frequency, frequency of maximum energy, duration and inter-pulse interval. It is commonly acknowledged that spectrograms do not give sufficient information to identify individual species within some genera (Russo and Jones, 2002). In these cases, we ascribed the calls to a group of species (this is the case for the genera Nyctalus/Eptesicus, Myotis and Plecotus). Therefore, we ran eight models at the species level and three at the genus level.

Modelling bat distribution
We recorded climate variables that could affect bat distribution in Chelsa V1.2 for the period 2006-2015 (Karger et al., 2017a, b). Since drought constrains primary productivity (and the concomitant availability of insects) in Mediterranean habitats (Nahal, 1981), we selected Mean Temperature of Warmest Quarter (bio10, TEMP) and Precipitation of Driest Month (PREC, bio14) as two potential drivers of bat distribution. In addition, since bats use open and tree-covered landscapes (Dietz and Kiefer, 2016), we included bare ground (BARE) and tree cover (TREE, including all tree species; percentage) from the Vegetation Continuous Fields MOD44B (GlobCover 2.2; Bicher on et al., 2008). Finally, we included elevation because it is a topographical feature with a different meaning than climatic variables. Although the elevation is strongly related to the climate in mountain areas (Pepin and Lundquist, 2008), we included this parameter as a complementary indicator of certain meteorological events (wind, sudden changes in temperature, storms) which can affect bat occurrence. Layers were used at a resolution of 30 arc sec.
The presence/absence of individual species and the explanatory variables and their squares (standardized to mean of zero and standard deviation of one) at the 274 sampling points were used to run the multGLM function (family "binomial") in the fuzzySim package (R Core Team 2017, version 3.4.1; Barbosa 2015). This function includes a stepwise forward-backward selection of variables and provides the best model according to Akaike's Information Criterion (AIC). It then produces favourability functions that remove the effect of species prevalence on the probabilities predicted by current logistic regressions . We used the Area Under Curve (AUC) for checking the models' performance (package 'modEvA'; Barbosa et al., 2016). AUC above 0.5 means that the model performs better than random and above 0.75 that has very good discrimination capacity (Fan et al., 2006). We then used the getPreds function (in fuzzySim package) to extrapolate the models to the whole study region (Fig. 1). Finally, we obtained a fuzzy richness value (richness herein) by adding up the favourability values of species and calculated a rarity index (rarity herein) following Estrada et al. (2011). This index increases as the occurrence probability of rare species increases. Cartographic data were managed with QGIS 2.18.11 (QGIS Development Team, 2016).

Gap analysis
The geographical patterning of the two indices were compared with the distribution of the protected areas reported by the Spanish Government (https://www.miteco.gob.es/es/biod iversidad/servicios/banco-datos-naturaleza/informacion-dispo nible/ENP.aspx, 2020). The protected areas were classified as National Park (Guadarrama National Park; NP), buffer area of the National Park (BNP), other protected reserves (Other; all of them natural parks and reserves and within the Natura 2000 network) and areas outside the network of protected areas (No protection). To assess the distribution of the study indices, we generated 985 randomly selected points within the study area that were distributed among areas with different conservation statuses. We used Kruskal-Wallis ANOVAs (classification factor: protection status) to test for differences in the study indices among protected and unprotected areas. The analyses were carried out with Rcmdr 3.4.1 (Fox and Bouchet-Valat, 2019).

Urban gradient
To investigate the potential effect of urbanization on the distribution of the study indices, we used the scores provided by the Human Footprint, an index of population density, human land use and infrastructure (Sanderson et al., 2002). Using the 985 random selected points within the study area, we tested whether richness and rarity indices were negatively correlated to urban gradients by using a non-parametric rank order Spearman correlation analysis. In addition, to explore whether at similar urban indexes the protected areas were located in the best sites for bats according to their richness and rarity, we used a non-parametric Quade's test (Quade, 1967). This test regresses the ranks of the dependent variable (here, richness or rarity) against ranks of the covariate (HFP), then it fits a one-way analysis of variance to the regression residuals, using the grouping variable (protected vs. non-protected) as the factor. A significant result between groups (e.g., higher values in protected areas) will indicate that the best sites for bats tend to be in protected areas along the urban gradient. Statistical analyses were carried out with Rcmdr 3.4.1 (Fox and Bouchet-Valat, 2019).

Modelling bat distribution
We recorded 11 bat species or groups of species whose occurrence data (Appendix S2) were used to produce distribution maps by using the favourability functions that were extrapolated to the whole study area (Table 1; Fig. 2). Except for Rhinolophus hipposideros and Tadarida teniotis models, whose AUCs were around 0.5 (random), models reflected a suitable discrimination capacity (Table 1; Fan et al., 2006). As a rule, all bats reported the highest occurrence probabilities at midelevation from which they decreased according to their particular preferences (see, e.g., Barbastella barbastella vs. Pipistrellus pygmaeus; Fig. 2). Rhinolophus hipposideros (the rarest species) did not respond to any variable in the model and Pipistrellus pipistrellus (the most common species) was detected in the full set of sampling points which made it impossible to model its occurrence probability (see Appendix S2). As a consequence, these species were not included in the multi-species indices (richness and rarity indices). These indices showed a similar pattern with the highest scores in the Guadarrama Mountains and the surrounding piedmont (Fig. 3). Interestingly, richness and rarity indices were positively correlated (Spearman correlation coefficient, r = 0.95, P < 0.001, n: 985) suggesting that the most suitable sites were also the best ones for the rarest species.

Gap analysis
The gap analysis reported overlap between the most important areas for bats and the network of protected areas (Figs. 3 and 4). As a result, the mean scores of species richness and rarity within all the protected areas were significantly higher than the mean scores detected outside the network (Kruskal-Wallis test; richness index, protected: 3.53 AE 0.03, not protected: 2.93 AE 0.03, v 2 1,984 = 123.75, P < 0.001; rarity index, protected: 0.98 AE 0.01, not protected: 0.84 AE 0.01, v 2 1,984 = 152.28, P < 0.001). Within the set of protected areas, the Guadarrama National Park (NP) combined with its buffer area (BNP) showed higher scores than the rest of the protected areas (Other)

Urban gradient
The two indices were negatively correlated to the human footprint, suggesting a negative effect of infrastructures and urban encroachment on the distribution of the best areas for bats (human footprint vs. richness index, Spearman r = À0.57, P < 0.001; human footprint vs. rarity index, Spearman r = À0.38, P < 0.001, n = 985; Fig. 5). Within this gradient, the protected areas tended to occur in the most suitable areas for bat richness (Quade's test, F 1,983 = 17.30, P < 0.001) and rarity (F 1,983 = 17.67, P < 0.001; Fig. 5).

Modelling bat distribution
Although individual bat species have idiosyncratic distributions due to their particular habitat preferences (Dietz and Kiefer, 2016), bats seem to show similar trends in European mountain ranges, since most of them are particularly frequent at mid-elevations (Jaberg and Guisan, 2001;Charbonnier et al., 2016). The results in this study support this view since most bat species showed the highest favourability at mid-elevations of the Guadarrama Mountains, an area covered by an interspersed distribution of woodlands, scrublands and pasturelands; suitable habitats for a broad range of bat species (Dietz and Kiefer, 2016; Figs. 1, 2; Table 1). The models also suggest that most species avoided the bare agricultural areas of the lowlands and the urban areas around Madrid city (Figs. 1, 2). However, it is interesting to note that these models also suggest that the reported distribution of bats around the mountains varies per species. For instance, distribution of the tree-dwelling Western Barbastelle Barbastella barbastellus reflected its status as a mountain species (Russo et al., 2005), with its most suitable areas in the mature native Scots pinewoods, a northern tree species that in this region reaches its southernmost range edge (Sinclair et al., 1999); the most suitable areas for the Mediterranean Kuhl's Pipistrelle Pipistrellus kuhlii (Sachanowicz et al., 2006) were located across the Mediterranean dry sclerophyllous woodlands of the piedmont; and the potential habitats for the ubiquitous Common Pipistrelle Pipistrellus pipistrellus (Davidson-Watts and Jones, 2006) expanded across the entire region (Table 1). The distribution pattern of the most suitable sites for bat species in the Guadarrama Mountains was reflected in the distribution of the two indices (Fig. 3). More explicitly, medium elevation ranges in the mountains emerged as the best places for bats according to the study indices. This pattern agrees with the usual patterns of species richness in dry mountains where productivity, a main driver of species richness (Cusens et al., 2012), increases between dry-low and cold-high elevations (Rahbek, 2005;McCain, 2009). A similar pattern of increased species richness at intermediate elevations of the Guadarrama Mountains has already been observed in other taxonomic groups (e.g., ants and birds; Flores et al., 2018) supporting the role of this mountain range as a regional hotspot of biodiversity.

Gap analysis
Networks of protected areas are often designed to conserve certain charismatic species or habitats and they frequently ignore the protection of more cryptic or less popular groups Intercept and estimates for the variables included in the final model for each species. AUC scores reflect model discriminative performance (scores above 0.5 mean that the model performs better than random). Abbreviations: BARE, bare ground cover; BARE 2 , quadratic bare ground cover; ELEV, elevation; ELEV 2 , quadratic elevation; n, number of presences for each species recorded in 274 sampling points; PREC, precipitation of driest month; PREC 2 , quadratic precipitation of driest month; TEMP, temperature of warmest quarter; TEMP 2 , quadratic temperature of warmest quarter; TREE, tree cover; TREE 2 , quadratic tree cover.
Animal Conservation 25 (2022) (Bosso et al., 2016). However, the gap analysis carried out in this study showed that the mean scores of the two indices were higher in the protected areas and that the protected area network in central Spain overlapped most of the best areas for bats (Fig. 3). These results support previous studies reporting that the network of protected areas in Europe tends to protect the best areas for some bat species (Kerbiriou et al., 2018). Interestingly, the Guadarrama National Park showed higher richness and rarity indices than the rest of the protected areas (Fig. 4), a good result in terms of conservation since it is within the highest protection levels of Spanish law (Spanish law 42/2007, del Patrimonio Natural y de la Biodiversidad).

Urban gradients
Increasing urbanization is usually related to a depletion of biodiversity (Pautasso, 2007;Goddard et al., 2010). However, this does not always occur because both biodiversity and human settlements respond positively to increasing levels of primary productivity (Chown et al., 2003). This means that human populated areas occur in areas of potential high biodiversity while the surrounding wild areas tend to be less productive and thus avoided by humans. In this context, it is interesting to note that the best areas for bats in central Spain (Fig. 3) were located outside the most urbanized areas (Fig. 1). This is probably an idiosyncratic trait of the study area where the urban encroachment of Madrid (an administrative city located ad hoc in the geographical centre of Spain) has expanded along a radial road and railway network and is not mediated by climate and agricultural suitability (Dowling, 2016). As a result, the distant wilderness areas of the Guadarrama Mountains have been hardly affected by urban encroachment. This process is depicted by the negative relationship between richness and rarity indices and urbanization intensity within the study area (Fig. 5). In addition, this trend showed that, regardless of the interspersion of urbanized versus bat-suitable areas, the protected areas tended to be located in the best areas for bats. This result supports two main ideas: Firstly, that despite the occurrence of bats in urban areas of the study area (Tena et al., 2020a), the group tends to disappear as urban encroachment increases (Avila- Flores and Fenton, 2005). Secondly, it corroborates the view that at local scales, the protected areas of central Spain protect the most adequate activity areas for bats.

Conclusions and conservation prospects
This work shows that the distribution of protected areas in central Spain overlaps with most of the best activity areas for bats as defined by species distribution models. This result is strongly related to some idiosyncratic features of the study area such as the outstanding conservation importance of    Sorte and Jetz, 2010). In fact, some invertebrates of the Guadarrama Mountains have already suffered an elevational shift accompanying recent climate change during the last decades (e.g., butterflies; Wilson et al., 2005Wilson et al., , 2007. These changes in the distribution of woodlands and some potential prey species could affect the distribution of the most suitable foraging areas for bats, which could expand to the highest elevations and retreat in the piedmont. These trends could also be affected by the ongoing process of forest regeneration and shrub expansion typical of central Spain and other Mediterranean areas (Kuemmerle et al., 2016;Morales-Molino et al., 2017), which have benefited populations of other forest species (e.g., birds; Teller ıa, 2019) but can eventually suffocate biodiversity through the pervasive effect of extreme shrub encroachment and tree densification (Regos et al., 2016;Tena et al., 2020b). However, since all of these potential changes will occur within a protected area network, it seems important to carry out conservation programs to anticipate and mitigate the potential effects of changes on bat assemblages. Unfortunately, this reserve network has not carried out effective monitoring of bat populations and their habitats are poorly known at regional (De Paz et al., 2015) and global scales (Mickleburgh et al., 2002) causing them to remain relatively invisible to legislators, conservationists and managers.
Finally, it is interesting to comment on the utility and limitations of bioacoustic sampling (poor discrimination capacity for Rhinolophus hipposideros and Tadarida teniotis models; Loiselle et al., 2003;Buckman-Sewald et al., 2014) as a way to generate data to model bat distribution without the use of sometimes skewed information provided by museum and/or bibliographic records. We consider that the most obvious usefulness of this approach is the quick mapping of bat habitat suitability within the study areas . Since this information refers to activity areas, it seems a promising complementary approach to the location and protection of bat roosts (Medell ın et al., 2017). However, more studies are required to advance the application of call recording and distribution modelling in the detection and monitoring of important areas for bats (Razgour et al., 2011(Razgour et al., , 2016. For instance, it could be interesting to improve the call libraries used to determine bat species in bioacoustic sampling since they still show significant gaps (Nyctalus, Eptesicus and Myotis; Russo and Jones, 2002). It is also important to improve the quality of species distribution models by establishing sound biological links between the environmental predictors and the targeted bat species and testing the predictive ability of the models (Fois et al., 2018). Distribution of (a) richness and (b) rarity indices for a bat assemblage in central Spain, sampled with acoustic detectors between 2014 and 2016, along human footprint gradient assessed at 985 random points within the study area. Richness index is the sum of the favourability values of all species and groups; rarity assigns the highest scores to areas where rare species have high occurrence probabilities. Distance weighted least squares lines show the patterns for protected (green circles) and non-protected (white squares) areas.