Modelling bat distributions and diversity in a mountain landscape using focal predictors in ensemble of small models

Bats are important components of mammalian biodiversity and strong bioindicators, but their fine‐scale distributions often remain less known than other taxa (e.g., plants, birds). Yet as highly mobile species with multiple needs in the landscape, bats impose serious modelling challenges, such as advanced use of neighbourhood analyses. The aims of this study were to test the use of a designed sampling of bats for biodiversity and conservation assessments, and to find appropriate modelling solutions for providing nature practitioners with reliable potential bat distribution maps in a mountain area of high conservation interest.


| INTRODUC TI ON
Pressures on ecosystems such as land use change, habitat fragmentation and climate change are increasing (e.g., Foley et al., 2005;Pachauri et al., 2014;Maxwell, Fuller, Brooks, & Watson, 2016). Therefore, understanding the interaction between species and their environment is paramount to identify environmental or climatic factors driving the distribution of species and to provide reliable information on current and future suitable habitat to allow efficient conservation management (Guisan et al., 2013).
Species distribution models (SDMs, based on the quantification of the envelope of suitable habitats for a species; also called ecological niche (ENMs) or habitat suitability (HSMs) models; see Guisan et al., 2013) are useful tools to reach this goal (Franklin, 2010;Guisan & Thuiller, 2005;Guisan, Thuiller, & Zimmermann, 2017;Guisan & Zimmermann, 2000;Peterson et al., 2011). These models rely on presence-absence or presence-only (occurrence) data and a set of environmental variables to quantify the species-environment relationship and approximate the realized environmental niche (or part of it) of a species in a given area. This niche can then be projected in geographic space to obtain potential distribution of the species, possibly at different times (e.g., past, future). Nowadays, SDMs are applied across a large array of taxa and ecological fields (e.g., assessment of the impact of land use, climate or other environmental changes on species' distribution; support of conservation planning and reserve selection or the identification of unsurveyed sites for rare or elusive species; Zanini, Pellet, & Schmidt, 2009;Broennimann et al., 2014;Maggini et al., 2014;Guisan et al., 2017).
Bats are the second most-diverse mammalian order with ca. 1,300 species worldwide (Fenton & Simmons, 2015) and often account for an important part of a region's vertebrate diversity. Additionally, insectivorous bats are also important bioindicators (Jones, Jacobs, Kunz, Willig, & Racey, 2009; due to their position in the food web making them particularly sensitive to habitat fragmentation and land use changes (Meyer, Struebig, & Willig, 2016;Mickleburgh, Hutson, & Racey, 2002;Voigt & Kingston, 2016). Despite their ecological success and potential as umbrella species (Lisón, Sánchez-Fernández, & Calvo, 2015;Roberge & Angelstam, 2004), bat ecology and distributions remain less known (compared to birds for instance), hampering monitoring and conservation efforts. This knowledge deficit is mostly due to bats' nocturnal and elusive lifestyle, which makes them difficult to observe and identify in the field. In recent years, the development and widespread adoption of acoustic monitoring techniques has greatly improved our ability to collect data about this order, despite the lack of a standardized technique (Darras, Pütz, Rembold, & Tscharntke, 2016;MacSwiney, Cristina, Clarke, & Racey, 2008;Stahlschmidt & Brühl, 2012) and difficulty of precise species-level identification (Barclay, 1999;Obrist, Boesch, & Flückiger, 2004;Russo & Voigt, 2016). Improvement in both sampling efficiency and modelling technique leads to an encouraging increase in studies using SDMs for bats (for a review, see Razgour, Rebelo, Febbraro, & Russo, 2016).
However, SDM studies on bats remain scarce and often address disparate questions or single species such as disentangling niches of cryptic species (e.g., Sattler, Bontadina, Hirzel, & Arlettaz, 2007;Rutishauser, Bontadina, Braunisch, Ashrafi, & Arlettaz, 2012;Santos et al., 2014;Smeraldo et al., 2018), modelling rare species of conservation concern (e.g., Razgour, Hanmer, & Jones, 2011;Santos, Rodrigues, Jones, & Rebelo, 2013;Silva, Vieira, Silva, & Cassia Faria, 2018) or predicting the influence of climate change (e.g., Rebelo, Tarroso, & Jones, 2010;Pio et al., 2014;Carstens, Morales, Field, & Pelletier, 2018). Applications of SDMs to multiple or all species in a region at fine spatial resolution are still relatively rare due to limited data availability (but see , although such effort would be needed to support conservation planning (Vincent, 2017). In recent years, a promising approach to deal with limited number of occurrences (i.e., rare or difficult to detect species) and large sets of environmental variables was developed (Breiner, Guisan, Bergamini, & Nobis, 2015;Breiner, Nobis, Bergamini, & Guisan, 2018) opening new modelling perspectives for such elusive groups (i.e., difficult to survey). This is especially useful when modelling bats as SDMs for these species should incorporate a large number of landscape elements, such as distance to resources, land structures and land cover variables (Rainho & Palmeirim, 2011), to capture the complexity of the ecological niche of these highly mobile species (Jaberg & Guisan, 2001).
Most published studies based on SDMs use gridded environmental maps (e.g., climate, soil, land use) to analyse and predict the potential distribution of species. This might be appropriate for sessile organisms (e.g., plants) but might be problematic if modelling highly mobile species such as bats or birds at high spatial resolutions, whereas for highly mobile species, the information on available resources within home/foraging range is often ecologically more important than the individual value of a single cell/location (e.g., 100 m × 100 m; Guisan & Thuiller, 2005). The key to integrate this information into SDMs might be focal variables representing an average/overview of the surrounding cells within a predefined focal window (Pellet, Guisan, & Perrin, 2004). However, even with focal variables remains the question of the optimal spatial scales (i.e., size of the focal window), depending on the location of accessible habitats (Eigenbrod, Hecnar, & Fahrig, 2008) and the foraging range and behaviour of the species Krauel, Ratcliffe, Westbrook, & McCracken, 2018;Pellet et al., 2004;Razgour et al., 2011;Villalobos-Chaves, Spínola-Parallada, Heer, Kalko, & Rodríguez-Herrera, 2017;Zanini et al., 2009). Jaberg and Guisan (2001) aggregated several land cover and landuse classes within large-enough pixel units of 2.5 km × 2.5 km to ensure capturing all important habitats and landscape components for bats. However, Bellamy, Scott, and Altringham (2013) showed that better answers to this scale question can be given by using a series of univariate models with different focal window sizes revealing highly variable and species-specific responses in bats.
In this study, we took the challenge to improve our capacity to gather distributional data for bat species in under-sampled regions, such as mountains, and use these data to build SDMs for use in improved conservation planning exercises within key conservation areas. We illustrate the approach in the Swiss Alps of Vaud, a key conservation area according to key non-governmental conservation agencies (Lassen & Savoia, 2005;WWF Vaud & Pro Natura Vaud, 2015), as well as a priority area for disciplinary and transdisciplinary research (see http://rechalp.unil.ch), where bats were not yet used in conservation planning (e.g., not included in Ramel, 2018). We collected field data from May to August by passive acoustic survey and mist-netting for bats species in the Western Swiss Alps. A multi-scale approach involving the test of nested focal circular windows of varying radius within univariate models was then used to select the best scale of influence for each land cover/use variable.
Finally, ensemble of small models (ESMs), an approach optimized to model the distributions of rare or under-sampled species (Breiner et al., 2015), were built at a resolution of 100 m × 100 m, using the mixed set of topo-climatic, distance and optimal focal land use/cover variables (i.e., based on the optimal radius per land use/cover class and per species). Our aims were more specifically to: (a) provide robustly designed data about as many of the bat species occurring in this area as possible, (b) test the use of distance and focal variables to investigate and better predict species-habitat relationships and (c) use these data and focal analyses to provide optimized spatially explicit potential distribution maps for use in conservation planning.

| Study area
The study area is located in the western Swiss Alps (Canton de Vaud, 46°10 to 46°30′N; 6°50′ to 7°10′E, Figure 1) and represents both a priority area for interdisciplinary research at the University of Lausanne (http://rechalp.unil.ch; Von Daniken, Guisan, & Lane, 2014) and a priority area for conservation at the European level (Lassen & Savoia, 2005). It covers ca. 700 km 2 , which encompass an elevation gradient ranging from 372 m to 3,210 m a.s.l. However, for this study with bats, only the potentially forested areas up to 2,000 m were considered for sampling and modelling, ignoring the 40 km 2 area above the treeline that are less occupied by bats , as only a few bat species are adapted to survive in high mountains (Alberdi, Aizpurua, Aihartza, & Garin, 2014;Le Roux et al., 2017;Weier, Linden, Gaigher, White, & Taylor, 2017). The study area and its vegetation are considerably influenced by the topography and human activity with the plain (Rhône valley) being densely populated and intensively farmed while the subalpine areas are essentially shaped by touristic activity and more extensive agricultural exploitations constituting a mosaic of meadows, pastures and forest patches (for more information about this area, see also http://rechalp.unil.ch; Guisan, 2005;Pellissier et al., 2012;Pradervand, 2015).

| Species data
So far, the bats in the study area were poorly documented except for five field surveys from the last national red list assessment (Bohnenstengel et al., 2014) and a long-term alpine bird-ringing station where bat catches occur regularly. New data about bat species were collected during a field survey from May to mid-August 2016. As ease of detection and identification can vary between species depending on the method of survey (Barnhart & Gillam, 2014;Flaquer, Torre, & Arrizabalaga, 2007), passive acoustic sampling was combined with capture using mist nets (Bohnenstengel et al., 2014;Stahlschmidt & Brühl, 2012). Additionally, the few historical data available for the area were incorporated. Species with less than 15 records in the area were not considered further, as they were not suitable for modelling. Furthermore, to avoid spatial autocorrelation between occurrences from different data sources and detectors (see below), a minimal distance of 1 km between occurrence points of each species was applied and other data within this radius were discarded. With a minimum distance of 1 km between occurrences, none of the bat species except the most frequently observed (Pipistrellus pipistrellus) showed significant spatial autocorrelation (Morans I > 0.05, Supporting Information Figure S1).

| Historical data
Historical data came from the national database (InfoFauna: Swiss Center for Faunal Cartography, http://www.cscf.ch). These data F I G U R E 1 Location of the sites sampled by passive acoustic or by mist nets capture from May to August 2016 were obtained with various technique, including passive and active acoustic, mist-netting and reporting of bats by the public. All reliable data with an accurate geolocation collected since 1991 between May and August were considered.

| Capture
Capturing by mist nets was done at 32 sites scattered across the study area with an average of 60 m of mist nets (2.4 m height, black M-14 14 mm monofilament from Ecotone©, Poland) set up in places favourable for bat passes until 4 hr after the sunset. Every bat was identified to the species level, and basic morphometry was taken (sex, age, forearm length, weight, reproductive condition). All the captures were done under licensing, and all bats were quickly released on site.

| Acoustic
Sixty-five sites were sampled acoustically following a regular sampling stratified by three altitudinal classes (Hirzel & Guisan, 2002).
On each site, two passive acoustic detectors (BatCorder ® 3.1, ecoObs©, Germany) were set up within a radius of 250 m from the sampling site and with a minimal distance of 100 m between them.
Whenever possible, the two detectors were placed in different habitat types within a site (e.g., open grassland/forest) to account for habitat heterogeneity and differing preferences among bat species.
If a species was recorded by both acoustic detectors placed per site one of the occurrences was discarded to prevent pseudo-replication.
Due to the rugged topography in the area making it difficult to access certain locations, some sampling points had to be slightly displaced.
Each site was recorded during 3 to 6 continuous nights, recording from one hour before sunset to one hour after sunrise. On each site, at least three nights offered suitable condition for bat activity based on temperature, precipitation and wind speed over the area.
Every file was automatically classified with SonoChiro ® (Biotope©, France) to allow an efficient selection of sequence to check per site.
Because low consensus has been reached, especially concerning the treatment of the acoustic data (Adams, Jantzen, Hamilton, & Fenton, 2012;Russo & Voigt, 2016), we let down the quantitative approach and limited ourselves to a qualitative list which could be checked manually. Therefore, for each detector per site, a call of each species was manually identified (Barataud, 2015) Rouse, Haas, Schell, & Deering, 1974) and the variance of canopy height [m]. We considered nine land use/cover classes and measured the proportion of each in concentric neighbouring windows of increasing radius of 0.1, 0.2, 0.3, 0.4, 0.5, 1, 1.5, 2, 2.5, 3, 4 and 5 km-to assess species-specific scales of influence of each variable. These variables were selected because they were not too highly correlated (spearman correlation <0.7; using the threshold recommended by Dormann et al., 2013) and have been shown to be highly influential for mobile species and bats in particular (Jaberg & Guisan, 2001). All environmental data were calculated at 100 m × 100 m resolution and additional details can be found in Table 1 and Supporting Information Appendix S1.

| Univariate models
For each species, the best window radius of the focal variables (the nine land use/cover classes and NDVI) was identified by using series of univariate models (Bellamy et al., 2013). Each of the 1,680 univariate models (14 species × 12 windows sizes × 10 focal variables) was run in R 3.3 (R Core Team, 2017) with the package biomod2 (Thuiller, Georges, Engler, & Breiner, 2016), using maxent 3.3.3 (Phillips, Anderson, & Schapire, 2006). For each model, we conducted 50 split-sample cross-validation runs (80% calibration/20% evaluation) using five different set of 5,000 pseudo-absences (selected with a minimum distance of 1 km from presences and down-weighted to equal the presence; i.e., prevalence of 0.5; Ferrier, Drielsma, Manion, & Watson, 2002). In a first step, the best spatial scale for each focal variable and species was selected based on the area under the curve (AUC) of a receiver operating-characteristic (ROC) plot (Fielding & Bell, 1997). In a second step, environmental variables were checked for multicollinearity and highly correlated variables (Spearman correlation > |0.7|; Dormann et al., 2013) were pruned by retaining the variable with the highest AUC (Bellamy et al., 2013).

| Ensemble of small models
All species distributions were predicted with an ensemble of small models approach (ESMs; Lomba et al., 2010;Breiner et al., 2015;Breiner et al., 2018) using the packages ecospat (Di Cola et al., 2017) and biomod2 (Thuiller et al., 2016) in R 3.3 (R Core Team, 2017). The ESM approach calibrates every possible bivariate combination of the environmental variables (i.e., n(n − 1)/2 bivariate models per species, with n being the number of variables) and then averages all those small bivariate models into an ensemble prediction. In this study, we only used generalized linear models (GLMs) for the calibration of bivariate models, as using different modelling technique within the ESM framework was shown not to improve performance significantly (Breiner et al., 2015(Breiner et al., , 2018. For each ESMs, we run 100 split-samples (80% calibration/20% evaluation) using ten different set of 10,000 pseudo-absences (selected with a minimum distance of 1 km from presences and down-weighted to equal the presence; i.e., prevalence of 0.5 ;Ferrier et al., 2002). Finally, an ensemble prediction was created by averaging the 100 runs weighted by AUC (Fielding & Bell, 1997). Weights of the bivariate models with an AUC lower than 0.5, hence worse than random (counter-predictions), were set to 0. The ensemble prediction was evaluated with the Continuous Boyce Index (Boyce, Vernier, Nielsen, & Schmiegelow, 2002;Hirzel, Lay, Helfer, Randin, & Guisan, 2006) Figure S2). Similar to standard model averaging or multi-model inference (Burnham & Anderson, 2002;Claeskens & Hjort, 2008), the contribution of each variable was assessed by calculating the mean weight of all bivariate models including the variable of interest, rescaled by the variable with the highest weight to a range between 0 and 1, as commonly done (but with Akaike's weights) in multi-model inference (e.g., Burnham & Anderson, 2002). These final ensemble models were then used to predict the potential distribution of our bat species in the study area up to 2,000 m (i.e., not projecting into the not sampled highest elevations).

| Reclassification of maps
To facilitate the use of the habitat suitability maps for practitioners and to account for the high mobility and varying home ranges of the bat species, the habitat suitability predictions were thresholded into binary maps, then smoothened using a focal window and finally categorized into different habitat classes. Binary presence-absence maps were generated by applying a threshold discarding the lowest 10% of training locations, hence encompassing 90% of training presences (minimal predicted area, MPA90; Engler, Guisan, & Rechsteiner, 2004). To account for the high mobility of bats, the binary distribution maps were then smoothed with a focal window of the mean home/foraging range of each species, giving for each target cell the proportion of presence cells within home range distance. The home range of each bat species was either provided by telemetric studies or expert-based (Arthur & Lemaire, 2009;Dietz, Nill, & Helversen, 2009). The use of binary predictions to build these smoothened distribution maps, hence the proportion of presence cells within home range distance instead of a mean habitat suitability within home range distance, was preferred because of the high mobility of the bats and good abilities to select optimal habitat and avoid unsuitable one in a heterogeneously suitable environment. The smoothened maps were reclassified into four classes of habitat suitability (HS) as follows: Unsuitable: Less than 5% of presences in home range distance; Marginal: Between 5% and 25% of presences; Suitable:

| Data collection
The field campaign did not target specific species and allowed us to obtain data for 18 different species of various ecology with 14 species having enough records for modelling (

| Univariate models-scale selection
The univariate models, aimed at selecting the best scale (i.e., focal window radius) per species for each land use/cover variable, revealed no clear trend (Figure 2). No best scale showed up among variables and species, suggesting scale dependences to be species-specific ( Figure 2). This univariate approach therefore proved particularly useful for selecting the optimal scales to include, for each variable, in the final ESM for each species. The full results of univariate models per species are provided in Supporting Information Figure S3.

| Ensemble of small models-species distribution
The ESM predictions for the 14 bat species were generally good,  Figures S4-S6).
In general, the four Euclidean distance variables had the highest contribution across species, especially the distances to roads, water and forest ( Figure 5). As expected in such mountain landscapes, slope and temperature were the most influential topo-climatic variables and outperformed all focal land use/cover variables (whatever the radius of the focal window). Among the land use/cover variables selected by the univariate models, the ones linked with forest (forest edge, deciduous forest, canopy structure) were more important than the others. However, depending on the species, other variables mostly related to human land use/cover (alpine pastures and cultivated land) had a high influence on the species distribution (Supporting Information Figure S7). The variables expressing productivity (i.e., primary production by NDVI at various scales) performed poorly for all species.
TA B L E 2 Species records in the study area after correction for autocorrelation (minimum distance 1 km) Note. Historical data were collected between May and August since 1991. Capture and passive acoustic were collected between May and August 2016. Total represents the number of occurrences used for modelling after removing duplicates. The national red list status (RLS) of each species is given as CR: critically endangered; EN: endangered; LC: least concern; NT: near threatened; VU: vulnerable.
The predicted species richness was higher in the lowlands, close to forest and along the major streams in the area ( The niche estimates, models and associated predictions derived from these various sampling and analytical steps bear the potential to provide important additional information to better understand the distribution and ecology of these highly mobile-and somewhat elusive-species as well as support future conservation planning exercises, as performed recently in Ramel (2018; including many taxa but not bats).
Our study illustrates how a carefully  were not detected sufficiently to fit models. Following one season of sampling, the results presented here thus change the status of the study area from a poorly surveyed to a well-surveyed area for bats in Switzerland, and our models and distribution maps will be able to serve as useful tools for conservation in future biodiversity assessments.
Our study further contributed to improve our knowledge and understanding of bat species ecology and distribution, especially compared to the two previous bat modelling studies in Switzerland (Jaberg & Guisan, 2001;Sattler et al., 2007), both conducted at the national scale and at much coarser resolution and lower spatial and environmental variability than the present study, which was conducted at high resolution and encompassed elevation from the lowlands (ca. 375 m) up to the highest elevations for bats (around 2,000 m). This study also innovated compared to previous studies by following a systematic screening and assessment of the most important scales of influence of the different land use/cover classes important for each species: The approach proposed here combines focal window analyses of various radii then tested within series of univariate models (here using GLMs), with-for each land use/cover class-the radii yielding the highest explained deviance being used in the final model. The big advantage of this approach is that it allows keeping a fine resolution (here 100 m) for environmental variables acting at very local scale (e.g., topography, micro-climate) while being able to account for neighbouring influence of the landscape up to several kilometres, which allows capturing the influence of multiple components of the landscape acting at larger scale that are key for bat presence (Bellamy et al., 2013). This was previously often accounted simply by using larger cell sizes (e.g., 2.5 km × 2.5 km in Jaberg & Guisan, 2001) that therefore include these many landscape components but at the cost of losing the fine resolution. The importance of this approach is especially shown by our finding that different radius are selected for different land use/cover classes and that these are species-specific, allowing each model to capture the unique combination of environmental variables important to define the species' environmental niche (see Guisan et al., 2017).
As expected, the projected species richness was highest at lower elevations associated with favourable climatic conditions for most species. In general, most bat species were highly associated with structural features of the landscape, such as roads,  (Bellamy et al., 2013;Willcox, Giuliano, Watine, Mills, & Andreu, 2017). Furthermore, when lit, roads also attract insects and therefore favour those species not disturbed by lights, such as Pipistrellus, Nyctalus and Eptesicus spp. (Rowse, Lewanzik, Stone, Harris, & Jones, 2016;Rydell, 2006;Spoelstra et al., 2017).
Additionally, the endangered species B. barbastellus was strongly associated with hedges and broadleaf forests confirming studies about habitat preference obtained from radio tracking individual females (Kühnert, Schönbächler, Arlettaz, & Christe, 2016;Zeale, Davidson-Watts, & Jones, 2012). Most bat species were associated with water ways, but this was especially prominent for P. pygmaeus known to be associated with lake shores and wide rivers in the warmer regions of Switzerland (Sattler et al., 2007). However, both model predictions and field records indicate that P. pygmaeus is expanding its range into valley's side slopes up to 900 m ( Figure 4).
In contrast to most species associated with warmer lowland conditions, E. nilssonii is limited to the higher parts of the area and avoids the lowlands (Figure 4). E. nilssonii seems to occur widely in the Prealps and was detected in every site above 900 m in eleva-  (Haupt, 2005;. While most bat species occupy most of their suitable climatic space, some species seem to fill only a small part of their climatic niche (Figure 4), such as Plecotus spp. that is most likely restricted due to habitat fragmentation making it especially vulnerable to land use change (Bohnenstengel et al., 2014 Different species seem to perform differently under anthropogenic pressure, with some species which could benefit and some which seem to be hampered by human's influence on the landscape .

| Study limitations and future perspectives
It is important to keep in mind the heterogeneity of the occurrence data used in this study. For all three data sources (passive acoustics, capture and historic data), it is impossible to determine whether the occurrence represents a breeding site, F I G U R E 5 Boxplot of the variable contribution of each environmental variable averaged across the 14 species. The percentage represents the proportion of the ensemble explained by the variable of interest. For variable explanation see Table 1.
The colours indicate the type of variable: red: Euclidean Distance; orange: topographic; blue: climatic; violet: focal; green: productivity F I G U R E 6 Species richness based on the sum of the binary predictions of each species (bS-SDM). On the left for all 14 species, in the middle for the nine species listed as LC or NT and on the right the five species listed as VU or EN in the national red list a foraging place or simply a migration corridor. We therefore can only model the probability to detect a certain species (Kery, Gardner, & Monnerat, 2010). We tried to account for this problem by smoothing our predictions with the individual home ranges of the different species, this way providing a smoother view of the species' range. However, future modelling efforts should preferably use occupancy models specifically developed to account for detectability (Kery et al., 2010), for instance using Bayesian approaches (Kery & Royle, 2015). As there is a large number of environmental and landscape factors potentially affecting the chance to detect a certain species, we decided to use ensemble of small models (ESM; Lomba et al., 2010;Breiner et al., 2015). The big advantage of ESMs is that even with a limited number of occurrences, a large number of predictors could be used (Breiner et al., 2015(Breiner et al., , 2018, but as to date, we do not know any development of occupancy modelling within an ESM context, which we see as a promising future development. Nevertheless, the models presented here still improved our understanding of the use of large and complex habitats by bats and could be used in future studies to predict the consequences of climate change and more direct anthropogenic alterations of the landscape, such as urbanization, roads construction, touristic activities, agricultural modifications or forest and water management, all included in a proper conservation planning prioritization scheme (see Tulloch et al., 2016). Additionally, our models and spatial predictions could help in identifying biodiversity hotspots, selecting the most valuable areas for protection and optimization conservation efforts of the most threatened places and species (e.g., Struebig, Christy, Pio, & Meijaard, 2010;Guisan et al., 2013;Pio et al., 2014), or be used to set up spatial monitoring programs for bats (Honrado, Pereira, & Guisan, 2016).

ACK N OWLED G EM ENTS
First and foremost, we like to thank David Progin whose Master thesis work provided most of the bat data and models used in this study, and contributed largely to the initial version of the manuscript. We