**2. Material and methods**

#### **2.1. Species data**

method, implemented as the free software 'Maxent' [11] is particularly popular. It attempts to find the distribution of maximum entropy (i.e., least constrained) that still agrees with all the observed data, and the value of the environmental variables at the locations where the species has been observed. Maxent performs well compared to other modeling methods [12], including when few presence data are available [13], making it especially attractive in data-poor regions. However, the method is vulnerable to bias in the input data [14]. It also shows a tendency to over fitting the presence data [15] and thus further enlarges the effect of sampling bias and spatial

In recent decades, a growing number of research studies have shown that niche models developed by incorporating remotely sensed predictors are more robust; these data can improve the prediction accuracy and tend to refine mapped distribution of species and habitats, com-

Remote sensing data can play an important role in developing cost-effective tools for modeling, mapping, planning, and protecting biodiversity. This is especially true at the scale of specific landscapes where the detection of patterns of species distribution can be greatly

In recent years, many new data on the distribution of bats in Bulgaria have been accumulated. Of particular importance in this respect was the project 'Mapping and Determining the Nature Conservation Status of Bats', activity 4, project DIR - 59,318-1-2 'Mapping and Determining the Nature Conservation Status of Natural Habitats and Species - Phase I', run between 2011 and 2013 by the Ministry of Environment and Waters. For complete project reports concerning bat species included in Annex 2 of the Habitat Directive see [18]. The abundant new data collected within the project, a result of intensive and extensive targeted studies for a brief period of time in context of the current state of nature, allow for a more in-depth analysis in the light of what has been known so far and the

This chapter aims to achieve the following: on the basis of presence-only modeling approach, combining current data on species distribution with a range of environmental layers, including satellite imagery, to reveal quantitatively the distributional patterns of bats in Bulgaria; to investigate potential ecological factors responsible for these patterns; to obtain meaningful biogeographical species groups; to document geographic patterns of species richness, rarity, and vulnerability; to analyze relationships between environmental factors, including anthropogenic changes of land cover and these biodiversity indices; and to highlight critical areas for

The result of the study can be useful for guiding further strategic conservation decisions, to assist the elaboration of management plans and to form a base for formulating restrictions and regimes to be included in future management plans of Natura 2000 sites, and to evaluate the impact of plans and projects on habitats and species listed in the Habitats

pared with climatic/topographical variables-only models [16].

improved by including this type of data [17].

existing knowledge gaps.

bat conservation.

Directive.

autocorrelation.

42 Bats

The data used in this study come from a database of georeferenced records developed as part of the project 'Mapping and Determining the Nature Conservation Status of Bats', activity 4, project DIR - 59,318-1-2 'Mapping and Determining the Nature Conservation Status of Natural Habitats and Species - Phase I', deposited in the Ministry of Environment and Waters and available under request. The data were collected during the period 2011–2012 within the Natura 2000 network. Bats were caught by hand from their roosts or by using mist-nets placed at entrances of caves, galleries, and at rivers and streams. All determinations were based on captured individuals, following the field guide of Dietz & Von Helversen [19]. Many of the captured individuals were photographed. Doubtful determinations were considered if their photographs and recorded standard measurements allowed the confirmation of the initial species identification. Exceptions were the determinations of the four species of the *M. mystacinus* morpho complex (*Myotis mystacinus/M. alcathoe* and *M. brandtii/M. aurascens*). Although the original determinations made by field experts were accepted, they should not be considered as certain, having in mind that 'further accumulation of genetic and morphologic data is needed to justify the variations and allow practical species identification' [20].

**Figure 1.** Spatial distribution of all georeferenced localities of 29 bat species used to model their habitat suitability. The number of points for each species is shown in **Table 1**.

In compiling the data set, the occurrence localities were screened to remove duplicate occurrences. In total, the final data set consisted of 1235 georeferenced point localities (**Figure 1**), which rendered 2766 unique species-locality combinations representing the distribution of 29 species. In Bulgaria, research efforts were mostly focused on bats roosting and swarming at caves and galleries. Relatively limited amount of records was available for the forest-dwelling species. The spatial coverage of the obtained data set clearly overcomes this discrepancy providing a balanced amount of data covering the bat's habitat diversity over the country.

Considering the importance of caves as roosts for a large part of the bats in Bulgaria, a digital layer was compiled, showing the presence of caves in a radius of 5 km (**Figure 2**). It was based

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived…

http://dx.doi.org/10.5772/intechopen.73623

45

Landsat imagery was used to present the peculiarities of the Earth's surface. In the selection of scenes, the aim was to represent the season of maximal vegetation development and to cover a period closest to the period during which the data were collected. However, in this regard, the author was struggling with some problems with the quality of the available images taken after 2003 due to a hardware failure and missing 22% of the pixels of Landsat 7 scenes [24]. For this reason, images taken before this date were selected. Ten scenes (georeferenced GeoTIFF files) from paths 181–185, rows 29–31 for the months of June–August of the years 1999–2000 from the Enhanced Thematic Mapper Plus (ETM+) sensor onboard Landsat 7 satellite were acquired from USGS [25] (for details see [26]). Despite the fact that the satellite images represent an earlier period, the peculiarities of the land cover showed the same general spatial patterns as during the study period. Correspondence between the state of the land cover presented by satellite imagery to that at the time of study is also because the data were collected within the Natura 2000 network, where the human impact on vegetation is poorly pronounced. This was also verified by comparing the land cover descriptions available in the

field protocols with that shown on the satellite images and on Google Earth™.

Adjacent bands in multispectral satellite images were highly correlated, which implies redundancy in the data. To overcome this problem, spectral values for different bands have been summarized by principal components analysis (PCA), and three uncorrelated principal components were extracted (**Figure 2**), presenting over 95% of the variability of the initial spectral data. To facilitate the interpretation of the environmental information presented by the obtained PCs, an analysis of the relationship between the principal component scores and the CORINE land cover classes at the third level [27] was made. For this purpose, for each cover type, the mean values of the pixels representing the principal component scores for each of the three components were calculated [25]. For the first Landsat PC1 (LPC1), the comparisons showed that moist places and forest vegetation (water bodies, watercourses, coniferous forest, mixed forest, broad-leaved forest) had the lowest pixel scores (average values of 1–1.5). Moors and heat land, transitional woodland shrub, green urban areas, sport and leisure facilities were associated with the middle part of the gradient (average values of 1.5–2). The range of 2–2.5 presented moderately anthropogenically influenced habitats such as agricultures and natural vegetation, fruit tree plantations, natural grassland, complex cultivation patterns, pastures, road and rail networks, vineyards, urban fabric, and nonirrigated arable land. The highest scores on this gradient (average values in the range of 2.5–3) represented xerophilous bare surfaces—sparsely vegetated areas, most of which are highly anthropogenically influenced—industrial units, mineral extraction sites, dump sites, dunes at beaches, and bare rocks. These comparisons show that LPC1 represents the nature of the vegetation—from welldeveloped forest vegetation in relatively humid conditions through shrubs and open spaces to bare and anthropogenically disturbed places. In practice, it can be said that it largely reflects the type of vegetation mainly in the form of the degree of its anthropogenic disturbance. Comparisons with the second principal component of spectral data (LPC2) showed that it separates the coniferous forests (low mean scores) compared to other types of vegetation, mainly to deciduous forests

on the coordinates of the studied caves, available in the database.

#### **2.2. Ecogeographical predictors**

Two topographical variables were used as proxies for abiotic conditions are as follows: elevation above sea level and topographic wetness index (**Figure 2**). The first variable was based on 1-arcsecond (30 m) SRTM digital elevation model freely available at the United States Geological Survey's EarthExplorer site [21]. Preliminary analyses have shown that on the surveyed territory, the elevation is strongly negatively correlated with the mean annual temperature, maximum temperature warmest month and quarter, minimum temperature coldest month and quarter, mean temperature wettest and driest quarter, while rainfall is positively correlated along the altitudinal gradient.

The digital elevation model was used to calculate topographic wetness index (TWI) using, SAGA GIS [22], which represents areas with high water retention potential (higher scores). It describes the depressions of the relief, where water bodies and river beds most often occur, and the local humidity is higher. It can be supposed that it will be an important ecogeographical predictor of habitat suitability for many bat species, in particular, in lower parts of the country which are hot and dry during summer. It is also known that many species of bats prefer the proximity of water and/or riverine vegetation where insect abundance tends to be higher [23].

**Figure 2.** Environmental layers. (A) Altitude. (B) Presence of caves. (C) Topographic wetness index. (D) First principal component of Landsat spectral data (LPC1). (E) Second principal component of Landsat spectral data (LPC2). (F) Third principal component of Landsat spectral data (LPC3).

Considering the importance of caves as roosts for a large part of the bats in Bulgaria, a digital layer was compiled, showing the presence of caves in a radius of 5 km (**Figure 2**). It was based on the coordinates of the studied caves, available in the database.

In compiling the data set, the occurrence localities were screened to remove duplicate occurrences. In total, the final data set consisted of 1235 georeferenced point localities (**Figure 1**), which rendered 2766 unique species-locality combinations representing the distribution of 29 species. In Bulgaria, research efforts were mostly focused on bats roosting and swarming at caves and galleries. Relatively limited amount of records was available for the forest-dwelling species. The spatial coverage of the obtained data set clearly overcomes this discrepancy providing a balanced amount of data covering the bat's habitat diversity over the country.

Two topographical variables were used as proxies for abiotic conditions are as follows: elevation above sea level and topographic wetness index (**Figure 2**). The first variable was based on 1-arcsecond (30 m) SRTM digital elevation model freely available at the United States Geological Survey's EarthExplorer site [21]. Preliminary analyses have shown that on the surveyed territory, the elevation is strongly negatively correlated with the mean annual temperature, maximum temperature warmest month and quarter, minimum temperature coldest month and quarter, mean temperature wettest and driest quarter, while rainfall is positively

The digital elevation model was used to calculate topographic wetness index (TWI) using, SAGA GIS [22], which represents areas with high water retention potential (higher scores). It describes the depressions of the relief, where water bodies and river beds most often occur, and the local humidity is higher. It can be supposed that it will be an important ecogeographical predictor of habitat suitability for many bat species, in particular, in lower parts of the country which are hot and dry during summer. It is also known that many species of bats prefer the proximity of water and/or riverine vegetation where insect abundance tends to be higher [23].

**Figure 2.** Environmental layers. (A) Altitude. (B) Presence of caves. (C) Topographic wetness index. (D) First principal component of Landsat spectral data (LPC1). (E) Second principal component of Landsat spectral data (LPC2). (F) Third

**2.2. Ecogeographical predictors**

44 Bats

correlated along the altitudinal gradient.

principal component of Landsat spectral data (LPC3).

Landsat imagery was used to present the peculiarities of the Earth's surface. In the selection of scenes, the aim was to represent the season of maximal vegetation development and to cover a period closest to the period during which the data were collected. However, in this regard, the author was struggling with some problems with the quality of the available images taken after 2003 due to a hardware failure and missing 22% of the pixels of Landsat 7 scenes [24]. For this reason, images taken before this date were selected. Ten scenes (georeferenced GeoTIFF files) from paths 181–185, rows 29–31 for the months of June–August of the years 1999–2000 from the Enhanced Thematic Mapper Plus (ETM+) sensor onboard Landsat 7 satellite were acquired from USGS [25] (for details see [26]). Despite the fact that the satellite images represent an earlier period, the peculiarities of the land cover showed the same general spatial patterns as during the study period. Correspondence between the state of the land cover presented by satellite imagery to that at the time of study is also because the data were collected within the Natura 2000 network, where the human impact on vegetation is poorly pronounced. This was also verified by comparing the land cover descriptions available in the field protocols with that shown on the satellite images and on Google Earth™.

Adjacent bands in multispectral satellite images were highly correlated, which implies redundancy in the data. To overcome this problem, spectral values for different bands have been summarized by principal components analysis (PCA), and three uncorrelated principal components were extracted (**Figure 2**), presenting over 95% of the variability of the initial spectral data. To facilitate the interpretation of the environmental information presented by the obtained PCs, an analysis of the relationship between the principal component scores and the CORINE land cover classes at the third level [27] was made. For this purpose, for each cover type, the mean values of the pixels representing the principal component scores for each of the three components were calculated [25]. For the first Landsat PC1 (LPC1), the comparisons showed that moist places and forest vegetation (water bodies, watercourses, coniferous forest, mixed forest, broad-leaved forest) had the lowest pixel scores (average values of 1–1.5). Moors and heat land, transitional woodland shrub, green urban areas, sport and leisure facilities were associated with the middle part of the gradient (average values of 1.5–2). The range of 2–2.5 presented moderately anthropogenically influenced habitats such as agricultures and natural vegetation, fruit tree plantations, natural grassland, complex cultivation patterns, pastures, road and rail networks, vineyards, urban fabric, and nonirrigated arable land. The highest scores on this gradient (average values in the range of 2.5–3) represented xerophilous bare surfaces—sparsely vegetated areas, most of which are highly anthropogenically influenced—industrial units, mineral extraction sites, dump sites, dunes at beaches, and bare rocks. These comparisons show that LPC1 represents the nature of the vegetation—from welldeveloped forest vegetation in relatively humid conditions through shrubs and open spaces to bare and anthropogenically disturbed places. In practice, it can be said that it largely reflects the type of vegetation mainly in the form of the degree of its anthropogenic disturbance. Comparisons with the second principal component of spectral data (LPC2) showed that it separates the coniferous forests (low mean scores) compared to other types of vegetation, mainly to deciduous forests (high mean scores). Open places and bushes occupy an intermediate position. In mountainous areas, it very well reflected the belts of forest vegetation. The range of scores on the third principal component (LPC3) was too narrow, which means that it reflected only slight differences in the nature of the land cover, which makes its interpretation difficult. Still, it can be said that it separated the sites covered with vegetation (lower mean scores) from those without vegetation, including most places of anthropogenic origin (higher mean scores). The environmental data set used in this study can be considered as sufficiently representative to characterize the ecological requirements of species. This provides a good basis for ecological modeling and it can be expected that the resulting models have the necessary ecological realism.

To summarize the results from modeling, a principal component analysis (PCA) was carried out on the matrix of correlations calculated from the 29 species models, using the 'Principal components' feature in ArcGIS v10.1 (ESRI). PCA summarizes the species models into a few new axes that best explain the variation found in the initial data set. The pixel scores of the first three axes, explaining the maximum variation were presented as raster layers. In order to figure out what these axes represent, the correlations between each of the PCA axes and the species raster files used as input were calculated. Similar correlations were also calculated between principal components and the environmental layers. Band Collection Statistics tool of ArcGIS v10.1 (ESRI) was used. The results were presented graphically in the form of ordination diagrams. Furthermore, based on the matrix of Euclidean distances between species in the coordinate system of the first three ordination axes, a dendrogram reflecting the similarity between the models was built. Ward's method was used for this

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived…

http://dx.doi.org/10.5772/intechopen.73623

47

Modeling data were further used to evaluate several synthetic parameters such as species richness, rarity, and vulnerability. For this purpose, the habitat suitability models were rescaled to a coarser resolution of 2.5 × 2.5 km. This grain can be considered close to the 'natural scale of resolution' for a species distribution model [33]. Various telemetry studies found

Two types of data were used for this purpose: quantitative, that is, the continuous Maxent habitat suitability values and qualitative, based on presence/absence values. In order to convert the continuous model predictions to discrete the presence/absence values, an arbitrary threshold of 50% was set. This threshold is too restrictive and probably lowered the actual number of species in a particular area. On the other hand, it corresponds to the recent recommendation to favor restrictive thresholds when stacking binary models to compute species richness [37]. Additionally, bearing in mind that the data were collected over all seasons, this threshold may help the regular presence of a species in an area (in summer, related to reproduction and in winter, related to hibernation) to be distinguished from accidental occurrence of vagrant or migrating individuals during spring and autumn [38]. Thus, it can be considered as representing the geographical locations of the excellent and optimal species habitats. Once the thresholds were set, the richness map was produced by combining/stacking binary

To develop the rarity pattern, we used the Index of Relative Rarity (IRR) normalized between 0 and 1 [39]. This index weighs species richness or total suitability according to the range sizes

maximum occurrences in the species pool, and r is the chosen rarity cut-off point (as percent-

× 0.97 + 1.05)

is the occurrence of species i, Qmin and Qmax are, respectively, the minimum and

2

(1)

or less [34–36].

purpose.

where Qi

age occurrence).

home ranges of species to encompass 2–3 km2

species maps, using the Raster Calculator feature in ArcGIS v10.1.

of the species present. Weights (W) were calculated as

<sup>W</sup> <sup>=</sup> exp(<sup>−</sup> <sup>Q</sup><sup>i</sup> \_\_\_\_\_\_\_\_\_\_\_ <sup>−</sup> <sup>Q</sup>min *ri* <sup>×</sup> <sup>Q</sup>max <sup>−</sup> <sup>Q</sup>min

Since the grain of Landsat and topography images (30 × 30 m) was too fine with respect to the scale at which the bat species discriminate and utilize essential habitat resources, environmental layers were rescaled to pixel size of 200 m. Although this resolution is still too fine, at least for climatic features of given area, it ensures a closer link between the values of other important factors and the bat locations. It corresponds well to such microhabitat features as the presence of suitable roost sites and foraging habitat [28].

#### **2.3. Statistical model**

In order to reveal the patterns of distribution of the species in the context of the available locations and the environmental factors considered to be relevant in this respect, modeling has been done by using the software Maxent (version 3.2.1). Model calculations were made using logistic output. Five cross-validation replicates were run for each species model and averaged into a single model. Recommended default values were used [11]. For each species, presence points were randomly divided into calibration (training) and evaluation (test) sets (30% samples for evaluation), and ROC curves and AUC figures were obtained. The potential sampling bias was addressed by the inclusion of the so-called bias file that allows generating background data with the same bias as occurrence data [29, 30]. The bias file was constructed using a Maxent model built on the basis of all locations. Model performance was evaluated based on the area under the curve (AUC) of the receiver operator characteristics (ROC) value. The AUC value is an indicator of the predictive accuracy of a model, correctly ranking presence locations higher than random. The AUC value ranges between 0 and 1, with higher values indicating better model fit; a model with an AUC = 0.5 indicates that the model performed no better than random, and a value over 0.75 is considered to be a good model performance [12]. The difference between the AUC values based on training and test localities was regarded as a measure of the degree of model over fitting. The smaller the difference between the two, the lesser the overfitting present in the model, resulting from preferential sampling and spatial autocorrelation [31].

#### **2.4. Postprocessing**

Based on obtained species models, Levin's niche breadths were calculated by using ENTools package. The measure treats suitability scores as proportional to utilization and reflects species' relative spread across a niche and vary in terms of suitability scores across space [32].

To summarize the results from modeling, a principal component analysis (PCA) was carried out on the matrix of correlations calculated from the 29 species models, using the 'Principal components' feature in ArcGIS v10.1 (ESRI). PCA summarizes the species models into a few new axes that best explain the variation found in the initial data set. The pixel scores of the first three axes, explaining the maximum variation were presented as raster layers. In order to figure out what these axes represent, the correlations between each of the PCA axes and the species raster files used as input were calculated. Similar correlations were also calculated between principal components and the environmental layers. Band Collection Statistics tool of ArcGIS v10.1 (ESRI) was used. The results were presented graphically in the form of ordination diagrams. Furthermore, based on the matrix of Euclidean distances between species in the coordinate system of the first three ordination axes, a dendrogram reflecting the similarity between the models was built. Ward's method was used for this purpose.

(high mean scores). Open places and bushes occupy an intermediate position. In mountainous areas, it very well reflected the belts of forest vegetation. The range of scores on the third principal component (LPC3) was too narrow, which means that it reflected only slight differences in the nature of the land cover, which makes its interpretation difficult. Still, it can be said that it separated the sites covered with vegetation (lower mean scores) from those without vegetation, including most places of anthropogenic origin (higher mean scores). The environmental data set used in this study can be considered as sufficiently representative to characterize the ecological requirements of species. This provides a good basis for ecological modeling and it can be expected

Since the grain of Landsat and topography images (30 × 30 m) was too fine with respect to the scale at which the bat species discriminate and utilize essential habitat resources, environmental layers were rescaled to pixel size of 200 m. Although this resolution is still too fine, at least for climatic features of given area, it ensures a closer link between the values of other important factors and the bat locations. It corresponds well to such microhabitat features as

In order to reveal the patterns of distribution of the species in the context of the available locations and the environmental factors considered to be relevant in this respect, modeling has been done by using the software Maxent (version 3.2.1). Model calculations were made using logistic output. Five cross-validation replicates were run for each species model and averaged into a single model. Recommended default values were used [11]. For each species, presence points were randomly divided into calibration (training) and evaluation (test) sets (30% samples for evaluation), and ROC curves and AUC figures were obtained. The potential sampling bias was addressed by the inclusion of the so-called bias file that allows generating background data with the same bias as occurrence data [29, 30]. The bias file was constructed using a Maxent model built on the basis of all locations. Model performance was evaluated based on the area under the curve (AUC) of the receiver operator characteristics (ROC) value. The AUC value is an indicator of the predictive accuracy of a model, correctly ranking presence locations higher than random. The AUC value ranges between 0 and 1, with higher values indicating better model fit; a model with an AUC = 0.5 indicates that the model performed no better than random, and a value over 0.75 is considered to be a good model performance [12]. The difference between the AUC values based on training and test localities was regarded as a measure of the degree of model over fitting. The smaller the difference between the two, the lesser the overfitting present in the model, resulting from preferential sampling and spatial

Based on obtained species models, Levin's niche breadths were calculated by using ENTools package. The measure treats suitability scores as proportional to utilization and reflects species' relative spread across a niche and vary in terms of suitability scores across space [32].

that the resulting models have the necessary ecological realism.

the presence of suitable roost sites and foraging habitat [28].

**2.3. Statistical model**

46 Bats

autocorrelation [31].

**2.4. Postprocessing**

Modeling data were further used to evaluate several synthetic parameters such as species richness, rarity, and vulnerability. For this purpose, the habitat suitability models were rescaled to a coarser resolution of 2.5 × 2.5 km. This grain can be considered close to the 'natural scale of resolution' for a species distribution model [33]. Various telemetry studies found home ranges of species to encompass 2–3 km2 or less [34–36].

Two types of data were used for this purpose: quantitative, that is, the continuous Maxent habitat suitability values and qualitative, based on presence/absence values. In order to convert the continuous model predictions to discrete the presence/absence values, an arbitrary threshold of 50% was set. This threshold is too restrictive and probably lowered the actual number of species in a particular area. On the other hand, it corresponds to the recent recommendation to favor restrictive thresholds when stacking binary models to compute species richness [37]. Additionally, bearing in mind that the data were collected over all seasons, this threshold may help the regular presence of a species in an area (in summer, related to reproduction and in winter, related to hibernation) to be distinguished from accidental occurrence of vagrant or migrating individuals during spring and autumn [38]. Thus, it can be considered as representing the geographical locations of the excellent and optimal species habitats. Once the thresholds were set, the richness map was produced by combining/stacking binary species maps, using the Raster Calculator feature in ArcGIS v10.1.

To develop the rarity pattern, we used the Index of Relative Rarity (IRR) normalized between 0 and 1 [39]. This index weighs species richness or total suitability according to the range sizes of the species present. Weights (W) were calculated as

$$W = \exp\left(-\frac{\mathbf{Q}\_i - \mathbf{Q}\_{\text{min}}}{r\_i \times \mathbf{Q}\_{\text{max}} - \mathbf{Q}\_{\text{min}}} \times 0.97 + 1.05\right)^2\tag{1}$$

where Qi is the occurrence of species i, Qmin and Qmax are, respectively, the minimum and maximum occurrences in the species pool, and r is the chosen rarity cut-off point (as percentage occurrence).

Occurrence-based index of relative rarity (OIRR) was calculated as

$$OIRR = \frac{\frac{\sum w\_i}{S} - w\_{\min}}{\frac{\upsilon}{w\_{\min}} - w\_{\min}} \tag{2}$$

**Species N-locs Niche** 

*Rhinolophus ferrumequinum*

*Rhinolophus hipposideros*

ecogeographic variables to the Maxent model.

**breath**

**Training AUC**

**Test AUC**

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived…

*Barbastella barbastellus* 68 0.74 0.71 0.66 35.99 16.87 1.76 21.78 23.11 0.48 *Eptesicus serotinus* 86 0.87 0.78 0.69 25.06 44.38 4.61 14.09 10.25 1.62 *Hypsugo savii* 87 0.84 0.75 0.70 55.37 13.90 10.22 6.80 7.05 6.67 *Miniopterus schreibersii* 145 0.69 0.80 0.71 15.62 63.36 11.36 0.82 5.44 3.41 *Myotis alcathoe* 30 0.42 0.79 0.71 56.94 2.98 13.09 19.20 6.80 0.99 *Myotis aurascens* 13 0.98 0.83 0.72 11.70 39.23 15.17 24.49 5.09 4.33 *Myotis bechsteinii* 65 0.89 0.80 0.73 6.75 14.01 11.73 22.75 32.59 12.17 *Myotis blythii* 108 0.90 0.79 0.74 34.87 36.80 10.68 5.00 8.45 4.20 *Myotis brandtii* 18 0.50 0.78 0.71 68.94 4.57 1.39 19.21 5.80 0.09 *Myotis capaccinii* 83 0.72 0.78 0.74 22.76 66.85 1.21 0.66 3.36 5.17 *Myotis daubentonii* 61 0.95 0.76 0.74 23.94 11.42 21.74 10.03 31.71 1.16 *Myotis emarginatus* 86 0.90 0.82 0.74 13.90 56.43 8.09 3.78 15.08 2.71 *Myotis myotis* 79 0.86 0.81 0.75 8.00 71.39 7.71 4.16 6.14 2.61 *Myotis mystacinus* 30 0.56 0.82 0.75 58.81 2.82 15.31 6.52 16.15 0.40 *Myotis nattereri* 50 0.80 0.85 0.77 12.92 58.98 1.74 8.41 9.55 8.41 *Nyctalus leisleri* 53 0.80 0.83 0.78 12.25 2.72 28.36 47.74 4.90 4.03 *Nyctalus noctula* 62 0.90 0.80 0.78 14.95 2.63 52.38 25.90 3.78 0.37 *Pipistrellus kuhlii* 6 0.70 0.85 0.79 72.84 17.58 1.82 0.00 5.50 2.27 *Pipistrellus nathusii* 6 0.89 0.88 0.79 11.96 17.22 40.12 1.80 0.07 12.17 *Pipistrellus pipistrellus* 72 0.95 0.86 0.82 20.20 1.48 22.11 34.09 19.70 2.42 *Pipistrellus pygmaeus* 17 0.92 0.87 0.82 0.00 4.98 2.13 65.63 18.97 8.29 *Plecotus auritus* 56 0.49 0.87 0.83 83.81 1.12 2.77 0.97 6.13 5.21 *Plecotus austriacus* 130 0.93 0.93 0.84 45.18 11.37 14.22 6.02 13.45 9.77 *Rhinolophus blasii* 28 0.61 0.87 0.85 44.52 36.87 12.98 0.18 3.48 1.98 *Rhinolophus euryale* 95 0.79 0.90 0.86 28.89 51.08 9.74 1.66 7.16 1.48

464 0.94 0.88 0.86 16.73 47.76 8.22 2.67 22.69 1.92

439 0.90 0.86 0.86 19.00 62.91 10.03 0.80 6.34 0.93

*Rhinolophus mehelyi* 12 0.67 0.90 0.87 48.83 30.86 5.66 8.33 4.94 1.39 *Vespertilio murinus* 32 0.64 0.88 0.88 63.81 1.52 11.25 1.04 16.13 6.24 Mean — 31.76 28.23 12.33 12.32 10.98 3.83

**Table 1.** Mean training and test AUC's for the replicate runs and estimates of the relative contributions of the

**Alt Cave TWI LPC2 LPC1 LPC3**

http://dx.doi.org/10.5772/intechopen.73623

49

where wi is the weight of the i-th species in the assemblage, S is the assemblage species richness, wmin and wmax are the minimum and maximum weight, respectively.

Abundance-based IRR (AIRR) was calculated as

$$AIRR \quad = \frac{\frac{\sum a\_i w\_i}{N} - w\_{\text{min}}}{w\_{\text{max}} - w\_{\text{min}}} \tag{3}$$

where *ai* and *wi* are, respectively, the habitat suitability and weight of the i-th species in the grid cell, N is total habitat suitability of all species for the cell, and *wmin* and *wmax* are the minimum and maximum weights, respectively. The package 'rarity' was used [40].

For calculation of vulnerability index (V), species were ranked according to the five threat categories defined by the International Union for Nature Conservation (IUCN 2017) as: (1) insufficiently known (data deficient), (2) least concern (3) near threatened, (4) vulnerable, and (5) endangered. The formula was

$$V = \frac{\sum\_{i=1}^{k} v\_{ii}}{S\_i} \tag{4}$$

where vn is vulnerability rank of species i and S<sup>r</sup> is the richness of cell r.

In order to investigate which environmental factors best explained these synthetic indices, multiple regression analysis was used. It can be expected that this approach violates the assumption that residuals should be independent and identically distributed, resulting in inflated type-I errors due to residual spatial autocorrelation [41]. However, recently, it has been shown that short-distance residual spatial autocorrelation and the associated inflated type-I errors have no significant influence on the interpretation of regression coefficients [42, 43]. Regression analyses were performed by using Statistica 7.0 software (StatSoft, Inc. Statistica for Windows, Tulsa, OK).

#### **3. Results**

#### **3.1. Distribution models**

All models had a high level of predictive accuracy (**Table 1**), with AUC-training values between 0.7 and 0.97. AUC-test values ranged between 0.65 and 0.9. The mean difference across species between AUC-training and AUC-test was 0.054 indicating some over fitting, that is, that some models fit too tightly to calibration data, limiting to some extent their ability to predict

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived… http://dx.doi.org/10.5772/intechopen.73623 49


Occurrence-based index of relative rarity (OIRR) was calculated as

ness, wmin and wmax are the minimum and maximum weight, respectively.

mum and maximum weights, respectively. The package 'rarity' was used [40].

\_\_\_\_ ∑*wi <sup>S</sup>* \_\_\_\_\_\_\_\_\_ <sup>−</sup> *<sup>w</sup>*min *w*max − *w*min

<sup>∑</sup>*<sup>a</sup>* \_\_\_\_\_ *<sup>i</sup> wi <sup>N</sup>* \_\_\_\_\_\_\_\_\_ <sup>−</sup> *<sup>w</sup>*min *w*max − *w*min

grid cell, N is total habitat suitability of all species for the cell, and *wmin* and *wmax* are the mini-

For calculation of vulnerability index (V), species were ranked according to the five threat categories defined by the International Union for Nature Conservation (IUCN 2017) as: (1) insufficiently known (data deficient), (2) least concern (3) near threatened, (4) vulnerable, and

> ∑ *i*=1 *S vri* \_\_\_\_\_ *Sr*

In order to investigate which environmental factors best explained these synthetic indices, multiple regression analysis was used. It can be expected that this approach violates the assumption that residuals should be independent and identically distributed, resulting in inflated type-I errors due to residual spatial autocorrelation [41]. However, recently, it has been shown that short-distance residual spatial autocorrelation and the associated inflated type-I errors have no significant influence on the interpretation of regression coefficients [42, 43]. Regression analyses were performed by using Statistica 7.0 software (StatSoft, Inc.

All models had a high level of predictive accuracy (**Table 1**), with AUC-training values between 0.7 and 0.97. AUC-test values ranged between 0.65 and 0.9. The mean difference across species between AUC-training and AUC-test was 0.054 indicating some over fitting, that is, that some models fit too tightly to calibration data, limiting to some extent their ability to predict

is the weight of the i-th species in the assemblage, S is the assemblage species rich-

are, respectively, the habitat suitability and weight of the i-th species in the

is the richness of cell r.

(2)

(3)

(4)

*OIRR* =

Abundance-based IRR (AIRR) was calculated as

*AIRR* =

*V* =

where vn is vulnerability rank of species i and S<sup>r</sup>

where wi

48 Bats

where *ai*

and *wi*

(5) endangered. The formula was

Statistica for Windows, Tulsa, OK).

**3. Results**

**3.1. Distribution models**

**Table 1.** Mean training and test AUC's for the replicate runs and estimates of the relative contributions of the ecogeographic variables to the Maxent model.

independent evaluation data. However, it should be mentioned that when the models were run without bias file (results not shown), the mean value of the differences in species models was 0.09.

**3.2. Species groups**

74.24% of the total variance.

degraded vegetation.

ecological groups (**Figure 5**).

Mediterranean species.

More information on the direction of impact of the ecogeographic variables used in modeling can be obtained from the results of the PCA based on the correlation matrix of the layers showing the habitat suitability of the species. The PCA correlation biplots showed distinct patterns in the distribution of bat species and ecogeographic variables (**Figure 3**). The first axis (spPC 1) was mainly described by altitude and explained 37% of total variation, whereas the second axis (spPC 2) was related to the presence of caves, contributing another 25.5% of total variation. The most important environmental descriptor for the third axis (spPC3) was LPC2 and to a lesser extent TWI. This axis explained 11.41% of the variability of the species models and represented the positive association of some species with the broad-leaved forests and the associations of some species with water bodies in lowlands. The three species axes explained

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived…

http://dx.doi.org/10.5772/intechopen.73623

51

The first principal component (spPC1), representing the influence of the altitude, and hence, the effect of the main climatic parameters, temperature and humidity, separated cold-lowing mountain species (high positive correlations) from some Mediterranean species occurring mainly in lowlands (high negative correlations with this axis) (**Figure 4**). The other species, albeit showing lower correlations with this axis, also demonstrated well-defined pattern with respect to altitude. Positive correlations on this axis exhibited species whose optima are in the mountain foothills and in the middle mountain belt (**Figures 3**, **4**). Poor correlations with this axis showed the species closely related to caves (**Figure 4**). However, their position on the ordination diagram indicated that they prefer the middle part of the elevation gradient (**Figure 3**). This is largely because most of the caves in Bulgaria are located in the foothills and lower parts of the mountains. Moderate negative correlations with the first axis showed species occurring mainly at the lower part of the elevation gradient (**Figures 3**, **4**). Regarding the second axis, cave-dwelling species formed a well-defined group in the uppermost part of the plot. Correlations to the third axis showed clearly the environmental preferences of species that are not closely confined to caves (**Figures 3**, **4**). They had positive correlations with this axis. These were the mesophilous species that form two groups on the ordination diagram (**Figure 3**)—those inhabiting the middle mountain belt, preferring deciduous forests (high positive correlations with LPC2) and those attached to the more open but humid habitats in the lower parts of the country. Species with negative correlations on this axis prefer dryer habitats with sparse vegetation. In this group, species associated with higher altitudes prefer rock cliffs, and those associated with the lower part of the elevation gradient prefer sites with

The dendrogram derived from the Euclidean distances between species within the coordinate system of the first three principal components showed the presence of several well-defined

Group 1 comprised relatively rare species with markedly discontinuous distribution pattern in Bulgaria occurring mainly in mountainous regions. With respect to habitat type, the majority of them are confined to moist and deciduous forests. *H. savii* is a petrophilous

Contrary to expectations, there was no correlation between the niche breadths and the AUCtest statistics. The same was valid for the number of points for each species. This contradicts to the existing opinions that the generalist species have lower AUC scores [12]. This result shows that even species with large niches, broadly distributed in the territory of the country are closely dependent on specific ecogeographic factors.

*Nyctalus noctula, Myotis emarginatus, M. blythii, Rhinolophus hipposideros, Pipistrellus pygmaeus, Plecotus austriacus, Rhinolophus ferrumequinum, Pipistrellus pipistrellus, Myotis daubentonii,* and *M. aurascens* had the widest niche. These are the species that are widespread throughout the country. With the narrowest niche were such species *as Myotis alcathoe, Plecotus auritus, Myotis brandtii, Myotis mystacinus*, confined to higher elevations.

The altitude and the presence of caves were the two ecogeographic variables that determine to the greatest extent the habitat suitability on the territory of the country. The altitude had the greatest impact on some mountain species, such as *Hypsugo savii, Myotis alcathoe, Myotis mystacinus, Vespertilio murinus, Myotis brandtii, Pipistrellus kuhlii, Plecotus auritus*, showing their strong commitment to the climatic peculiarities of the mountain regions in Bulgaria. The altitude hardly affected such species as *Pipistrellus pygmaeus, Myotis bechsteinii, Myotis myotis, Myotis aurascens, Pipistrellus nathusii, Nyctalus leisleri, Myotis nattereri, Myotis emarginatus, Nyctalus noctula*, which indicates that the climatic variables are of less significance compared to other factors.

The presence of caves determines to a large extent the suitability of the habitat and the distribution of some species such as *Eptesicus serotinus, Rhinolophus ferrumequinum, Rhinolophus euryale, Myotis emarginatus, Myotis nattereri, Rhinolophus hipposideros*, and *Myotis myotis*. The proximity of caves did not affect the habitat suitability of *Plecotus auritus, Pipistrellus pipistrellus, Vespertilio murinus, Nyctalus noctula, Nyctalus leisleri, Myotis mystacinus, Myotis alcathoe, Myotis brandtii, Pipistrellus pygmaeus*, and *Plecotus austriacus*.

The topographic wetness index reflecting the presence of water bodies in the lower parts of the country had the greatest impact on *Pipistrellus nathusii* and *Nyctalus noctula*.

The appearance of vegetation (coniferous vs. deciduous forests) as represented by LPC2 had a strong influence on dendrophilous species such as *Barbastella barbastellus, Myotis bechsteinii, Myotis aurascens, Nyctalus noctula, Pipistrellus pipistrellus, Nyctalus leisleri*, and *Pipistrellus pygmaeus*.

The effect of LPC1, representing the more general features of the land cover (wet places covered with predominantly dense forest vegetation vs. dry and bare sites) was less pronounced. This gradient had a strong influence on *Pipistrellus pygmaeus, Pipistrellus pipistrellus, Rhinolophus ferrumequinum, Barbastella barbastellus, Myotis daubentonii*, and *Myotis bechsteinii*. Most of them are dendrophilous species.

The effect of LPC3 was negligible—there were no species strongly affected by it.

#### **3.2. Species groups**

independent evaluation data. However, it should be mentioned that when the models were run without bias file (results not shown), the mean value of the differences in species models

Contrary to expectations, there was no correlation between the niche breadths and the AUCtest statistics. The same was valid for the number of points for each species. This contradicts to the existing opinions that the generalist species have lower AUC scores [12]. This result shows that even species with large niches, broadly distributed in the territory of the country

*Nyctalus noctula, Myotis emarginatus, M. blythii, Rhinolophus hipposideros, Pipistrellus pygmaeus, Plecotus austriacus, Rhinolophus ferrumequinum, Pipistrellus pipistrellus, Myotis daubentonii,* and *M. aurascens* had the widest niche. These are the species that are widespread throughout the country. With the narrowest niche were such species *as Myotis alcathoe, Plecotus auritus, Myotis* 

The altitude and the presence of caves were the two ecogeographic variables that determine to the greatest extent the habitat suitability on the territory of the country. The altitude had the greatest impact on some mountain species, such as *Hypsugo savii, Myotis alcathoe, Myotis mystacinus, Vespertilio murinus, Myotis brandtii, Pipistrellus kuhlii, Plecotus auritus*, showing their strong commitment to the climatic peculiarities of the mountain regions in Bulgaria. The altitude hardly affected such species as *Pipistrellus pygmaeus, Myotis bechsteinii, Myotis myotis, Myotis aurascens, Pipistrellus nathusii, Nyctalus leisleri, Myotis nattereri, Myotis emarginatus, Nyctalus noctula*, which indicates that the climatic variables are of less significance compared

The presence of caves determines to a large extent the suitability of the habitat and the distribution of some species such as *Eptesicus serotinus, Rhinolophus ferrumequinum, Rhinolophus euryale, Myotis emarginatus, Myotis nattereri, Rhinolophus hipposideros*, and *Myotis myotis*. The proximity of caves did not affect the habitat suitability of *Plecotus auritus, Pipistrellus pipistrellus, Vespertilio murinus, Nyctalus noctula, Nyctalus leisleri, Myotis mystacinus, Myotis alcathoe,* 

The topographic wetness index reflecting the presence of water bodies in the lower parts of

The appearance of vegetation (coniferous vs. deciduous forests) as represented by LPC2 had a strong influence on dendrophilous species such as *Barbastella barbastellus, Myotis bechsteinii, Myotis aurascens, Nyctalus noctula, Pipistrellus pipistrellus, Nyctalus leisleri*, and *Pipistrellus pygmaeus*.

The effect of LPC1, representing the more general features of the land cover (wet places covered with predominantly dense forest vegetation vs. dry and bare sites) was less pronounced. This gradient had a strong influence on *Pipistrellus pygmaeus, Pipistrellus pipistrellus, Rhinolophus ferrumequinum, Barbastella barbastellus, Myotis daubentonii*, and *Myotis bechsteinii*.

the country had the greatest impact on *Pipistrellus nathusii* and *Nyctalus noctula*.

The effect of LPC3 was negligible—there were no species strongly affected by it.

are closely dependent on specific ecogeographic factors.

*brandtii, Myotis mystacinus*, confined to higher elevations.

*Myotis brandtii, Pipistrellus pygmaeus*, and *Plecotus austriacus*.

Most of them are dendrophilous species.

was 0.09.

50 Bats

to other factors.

More information on the direction of impact of the ecogeographic variables used in modeling can be obtained from the results of the PCA based on the correlation matrix of the layers showing the habitat suitability of the species. The PCA correlation biplots showed distinct patterns in the distribution of bat species and ecogeographic variables (**Figure 3**). The first axis (spPC 1) was mainly described by altitude and explained 37% of total variation, whereas the second axis (spPC 2) was related to the presence of caves, contributing another 25.5% of total variation. The most important environmental descriptor for the third axis (spPC3) was LPC2 and to a lesser extent TWI. This axis explained 11.41% of the variability of the species models and represented the positive association of some species with the broad-leaved forests and the associations of some species with water bodies in lowlands. The three species axes explained 74.24% of the total variance.

The first principal component (spPC1), representing the influence of the altitude, and hence, the effect of the main climatic parameters, temperature and humidity, separated cold-lowing mountain species (high positive correlations) from some Mediterranean species occurring mainly in lowlands (high negative correlations with this axis) (**Figure 4**). The other species, albeit showing lower correlations with this axis, also demonstrated well-defined pattern with respect to altitude. Positive correlations on this axis exhibited species whose optima are in the mountain foothills and in the middle mountain belt (**Figures 3**, **4**). Poor correlations with this axis showed the species closely related to caves (**Figure 4**). However, their position on the ordination diagram indicated that they prefer the middle part of the elevation gradient (**Figure 3**). This is largely because most of the caves in Bulgaria are located in the foothills and lower parts of the mountains. Moderate negative correlations with the first axis showed species occurring mainly at the lower part of the elevation gradient (**Figures 3**, **4**). Regarding the second axis, cave-dwelling species formed a well-defined group in the uppermost part of the plot. Correlations to the third axis showed clearly the environmental preferences of species that are not closely confined to caves (**Figures 3**, **4**). They had positive correlations with this axis. These were the mesophilous species that form two groups on the ordination diagram (**Figure 3**)—those inhabiting the middle mountain belt, preferring deciduous forests (high positive correlations with LPC2) and those attached to the more open but humid habitats in the lower parts of the country. Species with negative correlations on this axis prefer dryer habitats with sparse vegetation. In this group, species associated with higher altitudes prefer rock cliffs, and those associated with the lower part of the elevation gradient prefer sites with degraded vegetation.

The dendrogram derived from the Euclidean distances between species within the coordinate system of the first three principal components showed the presence of several well-defined ecological groups (**Figure 5**).

Group 1 comprised relatively rare species with markedly discontinuous distribution pattern in Bulgaria occurring mainly in mountainous regions. With respect to habitat type, the majority of them are confined to moist and deciduous forests. *H. savii* is a petrophilous Mediterranean species.

**Figure 3.** Graphical representation of Pearson's correlation coefficients for the 29 bat species (circles) and ecogeographic variables (arrows) with the first three species PCA axes. (A) PC1 vs PC2. (B) PC1 vs PC3. The size of circles corresponds to the species niche breadth. ALT – Altitude, caves – Presence of caves within radius of 5 km, TWI – Topographic wetness index, LPC1–LPC3 – Landsat principal components. The species names are represented by three-letter abbreviations of the generic and species names.

Group 2. Two subgroups can be distinguished. Group 2a consisted of species with wide niches occurring at low, medium, and medium-high altitudes; although preferring caves, they also inhabit non-karstic areas where they use a variety of other roost types such as mine galleries, buildings, rock crevices, and hibernating in underground spaces. Group 2b comprised species more closely confined to caves, preferring karstic territories at the lower range of the altitudinal gradient.

**Figure 4.** Results of principal component analysis of species models. (A) Principal component scores on the first PC (sppPC1)—High scores suggest the habitat suitability for higher altitude species; the low ones represent the habitat suitability for warm-loving species in the lower parts of the country. (B) Principal component scores on the second PC (sppPC2)—High scores represent distribution of cave-dwelling bats. (C) Principal component scores on the third

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived…

http://dx.doi.org/10.5772/intechopen.73623

53

PC (sppPC3)—High scores suggest high habitat suitability for forest mesophilous species.

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived… http://dx.doi.org/10.5772/intechopen.73623 53

**Figure 4.** Results of principal component analysis of species models. (A) Principal component scores on the first PC (sppPC1)—High scores suggest the habitat suitability for higher altitude species; the low ones represent the habitat suitability for warm-loving species in the lower parts of the country. (B) Principal component scores on the second PC (sppPC2)—High scores represent distribution of cave-dwelling bats. (C) Principal component scores on the third PC (sppPC3)—High scores suggest high habitat suitability for forest mesophilous species.

Group 2. Two subgroups can be distinguished. Group 2a consisted of species with wide niches occurring at low, medium, and medium-high altitudes; although preferring caves, they also inhabit non-karstic areas where they use a variety of other roost types such as mine galleries, buildings, rock crevices, and hibernating in underground spaces. Group 2b comprised species more closely confined to caves, preferring karstic territories at the lower range of the altitudi-

**Figure 3.** Graphical representation of Pearson's correlation coefficients for the 29 bat species (circles) and ecogeographic variables (arrows) with the first three species PCA axes. (A) PC1 vs PC2. (B) PC1 vs PC3. The size of circles corresponds to the species niche breadth. ALT – Altitude, caves – Presence of caves within radius of 5 km, TWI – Topographic wetness index, LPC1–LPC3 – Landsat principal components. The species names are represented by three-letter abbreviations of

nal gradient.

52 Bats

the generic and species names.

**Figure 5.** Dendrogram showing the similarity between bat species based on their correlations with the first three principal components.

Group 3. Two subgroups were evident here. Group 3a comprised dendrophilous species, inhabiting wooded and humid areas in low and medium elevation. Group 3b embraced species restricted to lower altitudes, often associated with open areas, xerophilous sparsely vegetated rocky places or water basins.

**Figure 6.** Geographical patterns of bat species richness (A), rarity (B) presence/absence data, (C)-quantitative data, and

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived…

http://dx.doi.org/10.5772/intechopen.73623

55

**Table 2.** Standardized regression coefficient of multiple regression showing the relative contribution of each environmental variable in the prediction of species richness (SR), rarity (RA), and vulnerability (VU). All models are highly

**Environmental predictors SR Ra Vu** Cave PA 0.69 0.38 0.39 WI −0.12 0.17 −0.449 lpc3 0.06 0.02 −0.299 lpc2 0.22 −0.09 0.079 lpc1 −0.35 −0.18 −0.10 Alt 0.018 0.56 −0.40 Adjusted R2 0.82 0.59 0.46

vulnerability (D).

significant at P < 0.0001.

#### **3.3. Species richness, rarity, and vulnerability**

Superimposing the Maxent species models and derived weighted values according to rarity and vulnerability, resulted in the species maps of richness, rarity, and vulnerability, presented in **Figure 6**.

Correlation coefficients between these indices showed that the number of species was moderately positively correlated with rarity (0.67) and relatively less correlated with vulnerability (0.43), while the last two metrics were not correlated with each other. The beta coefficients of individual environmental factors as well as their overall contribution in the prediction of vulnerability (VU), species richness (SR), and rarity (RA) are presented in **Table 2**. The species richness was strongly affected by the presence of caves and moderately negatively influenced by LPC1. In total, the model explained 82% of the variance in species richness. Rarity, presenting the concentration of rare species was moderately positively affected by the altitude and to some extent by the presence of caves. Vulnerability was positively affected by presence of caves and negatively by TWI, altitude, and LPC3.

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived… http://dx.doi.org/10.5772/intechopen.73623 55

**Figure 6.** Geographical patterns of bat species richness (A), rarity (B) presence/absence data, (C)-quantitative data, and vulnerability (D).

Group 3. Two subgroups were evident here. Group 3a comprised dendrophilous species, inhabiting wooded and humid areas in low and medium elevation. Group 3b embraced species restricted to lower altitudes, often associated with open areas, xerophilous sparsely veg-

**Figure 5.** Dendrogram showing the similarity between bat species based on their correlations with the first three principal

Superimposing the Maxent species models and derived weighted values according to rarity and vulnerability, resulted in the species maps of richness, rarity, and vulnerability, pre-

Correlation coefficients between these indices showed that the number of species was moderately positively correlated with rarity (0.67) and relatively less correlated with vulnerability (0.43), while the last two metrics were not correlated with each other. The beta coefficients of individual environmental factors as well as their overall contribution in the prediction of vulnerability (VU), species richness (SR), and rarity (RA) are presented in **Table 2**. The species richness was strongly affected by the presence of caves and moderately negatively influenced by LPC1. In total, the model explained 82% of the variance in species richness. Rarity, presenting the concentration of rare species was moderately positively affected by the altitude and to some extent by the presence of caves. Vulnerability was positively affected by presence of

etated rocky places or water basins.

sented in **Figure 6**.

components.

54 Bats

**3.3. Species richness, rarity, and vulnerability**

caves and negatively by TWI, altitude, and LPC3.


**Table 2.** Standardized regression coefficient of multiple regression showing the relative contribution of each environmental variable in the prediction of species richness (SR), rarity (RA), and vulnerability (VU). All models are highly significant at P < 0.0001.

The results obtained from the modeling show that they correspond well to the knowledge about the ecology and the distribution of the species. The obtained result agrees with the prevailing opinions that, compared with less mobile species, the realized distributions of bats correspond closely with their potential distributions [44, 45]. This makes it clear that this approach is reliable and useful in two respects. On the one hand, it allows to identify the influence of individual environmental factors and to quantify their relative effect on species. Thus, modeling is also a realistic method for studying the ecology of individual species. On the other hand, on the basis of the statistical relations, it allows a reliable forecasting of species distribution.

areas of the country. The modeling results also show that in the lower parts of the country, the few available favorable habitats are too fragmented. This confirms the need to protect the

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived…

http://dx.doi.org/10.5772/intechopen.73623

57

The correlation of environmental predictors with species richness indicates that it is greatest in karstic areas and decreases in the direction of xerophilous, bare, and anthropogenically disturbed habitats, which is not a surprise. It is interesting to note, however, that species richness is not influenced by altitude (**Table 2**). It is often claimed that elevation gradient mirrors the latitudinal one, and species richness is assumed to decrease monotonically because of reduced temperature and consequent decrease in productivity. This is an indisputable fact on a larger geographic scale across continents. In fact, as it can be seen on the map (**Figure 6**), species richness appears to be higher in the middle mountain range. This is to a certain extent contradicts to the existing view that the species richness in Bulgaria is highest at the lower part of the elevational gradient, between 100 and 300 m, reaching 17–20 species [47]. This pattern can be explained by the fact that the role of changes in factors other than temperature is also significant. In the lower parts of the country, the forest is heavily destroyed, and the habitats of the bats are disturbed on large areas. Conversely, in the middle mountain belt, because of the rugged terrain, the anthropogenic disturbance of natural vegetation is weaker. Furthermore, elevation gradients have a more or less stable condensation zone at middle altitudes causing favorable conditions for many taxa, including invertebrates during the entire vegetation season (summer drought act in an opposite direction in lower elevations). On the variable terrain, local climate can vary considerably over short distances allowing small areas to present climatic optima for many bat species. Thus, middle mountain belt, being a transient zone, allows the co-occurrence of cold-loving species, preferring mountainous areas and some southern species, with wider ecological tolerance. In the highest parts of the mountains, species richness really decline. Thus, the change of the species richness of bats with altitude in this country is hump-shaped with a maximum at the lower part of the altitudinal gradient. This is one of the four common patterns noted in the literature, resulting from optimal combination between water availability and temperature [48, 49]. For bats in Bulgaria, however, the anthropogenic factor is also likely to be of significant importance for

Rarity primarily provides an insight into the facet of species biodiversity that is most at risk of extinction [50], also with respect to the maintenance of vulnerable ecosystem functions [51]. Different axes of rarity are usually considered: restricted abundance, restricted geographic distribution, and narrow niche breadth. In this study, the metric mainly reflects the spatial aspect. The obtained results suggest that, on the territory of Bulgaria, rare species occur predominantly in the karstic and mountain regions. In most cases, these are obviously species with limited tolerance with respect to the environmental conditions in the country. Of particular interest are the rare species confined to the mountains. It can be supposed that they will be most vulnerable to the climate change. It can be expected that the ecological amplitude for these species will continue to shrink, and their abundance will decrease during climate warming. Rare species associated with caves will obviously be vulnerable to human disturbances in caves and their

connecting landscape elements in the lower part of the country.

the observed trend.

surroundings.

The obtained results show that species form clusters that can be clearly explained in terms of analyzed ecogeographic variables. Generally, altitude is the key factor responsible for the largest differentiation between the species, especially between species in group 1 versus those in group 3b. The species placed in the center of PCA plots, belonging to group 2, tend to occur in environments characterized by the intermediate values of the investigated factors and hence explaining their wide distribution. The influence of altitude is not unexpected given that, in Bulgaria, it determines the spatial differentiation of the main climatic parameters such as temperature and rainfall [6]. Both parameters are known to impose limits on the ability of bats to forage for food [44, 45] and to survive prolonged hibernation periods or seasonal heat [46]. These climatic factors may also act indirectly by limiting the availability of essential resources.

Two factors related to geomorphology (presence of caves and topographic wetness index) are also considered to be important factors shaping bats distribution patterns. The presence of caves is important for the majority of the species. Although the topographical wetness index does not have a strong influence on most species, it largely determines the distribution of *Pipistrellus nathusii, Myotis daubentonii*, and *Nyctalus noctula* (**Figure 3**). The obtained result agrees well with the available knowledge that these species prefer the proximity of water and/ or riverine vegetation where insect abundance tends to be higher [4, 5].

Surprisingly, the effect of vegetation-related variables, LPC1 and LPC2, was relatively weak. Nevertheless, it should be mentioned that, LPC2, presenting the influence of deciduous forests, appears to be an important factor determining the habitat suitability of a large number of species, especially those that are not closely confined to caves. The relatively weak influence of this factor is largely a consequence of its correlation with the altitude and indirectly with the degree of anthropogenic disturbance of the forest cover, having in mind that most of the forests in Bulgaria are preserved mainly in mountains. Regardless of this correlation, the independent influence of this factor is very well demonstrated with respect to some species and allows outlining the distribution of some poorly known dendrophilous species. For some rare species, such as *M. bechsteinii* and *Nyctalus leisleri*, the positive influence of deciduous forests is much more pronounced than that of altitude. They also have large niche breadths, indicating that the rarity of these species is mainly due to the destruction of forests on large areas of the country. The modeling results also show that in the lower parts of the country, the few available favorable habitats are too fragmented. This confirms the need to protect the connecting landscape elements in the lower part of the country.

**4. Discussion**

56 Bats

distribution.

resources.

The results obtained from the modeling show that they correspond well to the knowledge about the ecology and the distribution of the species. The obtained result agrees with the prevailing opinions that, compared with less mobile species, the realized distributions of bats correspond closely with their potential distributions [44, 45]. This makes it clear that this approach is reliable and useful in two respects. On the one hand, it allows to identify the influence of individual environmental factors and to quantify their relative effect on species. Thus, modeling is also a realistic method for studying the ecology of individual species. On the other hand, on the basis of the statistical relations, it allows a reliable forecasting of species

The obtained results show that species form clusters that can be clearly explained in terms of analyzed ecogeographic variables. Generally, altitude is the key factor responsible for the largest differentiation between the species, especially between species in group 1 versus those in group 3b. The species placed in the center of PCA plots, belonging to group 2, tend to occur in environments characterized by the intermediate values of the investigated factors and hence explaining their wide distribution. The influence of altitude is not unexpected given that, in Bulgaria, it determines the spatial differentiation of the main climatic parameters such as temperature and rainfall [6]. Both parameters are known to impose limits on the ability of bats to forage for food [44, 45] and to survive prolonged hibernation periods or seasonal heat [46]. These climatic factors may also act indirectly by limiting the availability of essential

Two factors related to geomorphology (presence of caves and topographic wetness index) are also considered to be important factors shaping bats distribution patterns. The presence of caves is important for the majority of the species. Although the topographical wetness index does not have a strong influence on most species, it largely determines the distribution of *Pipistrellus nathusii, Myotis daubentonii*, and *Nyctalus noctula* (**Figure 3**). The obtained result agrees well with the available knowledge that these species prefer the proximity of water and/

Surprisingly, the effect of vegetation-related variables, LPC1 and LPC2, was relatively weak. Nevertheless, it should be mentioned that, LPC2, presenting the influence of deciduous forests, appears to be an important factor determining the habitat suitability of a large number of species, especially those that are not closely confined to caves. The relatively weak influence of this factor is largely a consequence of its correlation with the altitude and indirectly with the degree of anthropogenic disturbance of the forest cover, having in mind that most of the forests in Bulgaria are preserved mainly in mountains. Regardless of this correlation, the independent influence of this factor is very well demonstrated with respect to some species and allows outlining the distribution of some poorly known dendrophilous species. For some rare species, such as *M. bechsteinii* and *Nyctalus leisleri*, the positive influence of deciduous forests is much more pronounced than that of altitude. They also have large niche breadths, indicating that the rarity of these species is mainly due to the destruction of forests on large

or riverine vegetation where insect abundance tends to be higher [4, 5].

The correlation of environmental predictors with species richness indicates that it is greatest in karstic areas and decreases in the direction of xerophilous, bare, and anthropogenically disturbed habitats, which is not a surprise. It is interesting to note, however, that species richness is not influenced by altitude (**Table 2**). It is often claimed that elevation gradient mirrors the latitudinal one, and species richness is assumed to decrease monotonically because of reduced temperature and consequent decrease in productivity. This is an indisputable fact on a larger geographic scale across continents. In fact, as it can be seen on the map (**Figure 6**), species richness appears to be higher in the middle mountain range. This is to a certain extent contradicts to the existing view that the species richness in Bulgaria is highest at the lower part of the elevational gradient, between 100 and 300 m, reaching 17–20 species [47]. This pattern can be explained by the fact that the role of changes in factors other than temperature is also significant. In the lower parts of the country, the forest is heavily destroyed, and the habitats of the bats are disturbed on large areas. Conversely, in the middle mountain belt, because of the rugged terrain, the anthropogenic disturbance of natural vegetation is weaker. Furthermore, elevation gradients have a more or less stable condensation zone at middle altitudes causing favorable conditions for many taxa, including invertebrates during the entire vegetation season (summer drought act in an opposite direction in lower elevations). On the variable terrain, local climate can vary considerably over short distances allowing small areas to present climatic optima for many bat species. Thus, middle mountain belt, being a transient zone, allows the co-occurrence of cold-loving species, preferring mountainous areas and some southern species, with wider ecological tolerance. In the highest parts of the mountains, species richness really decline. Thus, the change of the species richness of bats with altitude in this country is hump-shaped with a maximum at the lower part of the altitudinal gradient. This is one of the four common patterns noted in the literature, resulting from optimal combination between water availability and temperature [48, 49]. For bats in Bulgaria, however, the anthropogenic factor is also likely to be of significant importance for the observed trend.

Rarity primarily provides an insight into the facet of species biodiversity that is most at risk of extinction [50], also with respect to the maintenance of vulnerable ecosystem functions [51]. Different axes of rarity are usually considered: restricted abundance, restricted geographic distribution, and narrow niche breadth. In this study, the metric mainly reflects the spatial aspect. The obtained results suggest that, on the territory of Bulgaria, rare species occur predominantly in the karstic and mountain regions. In most cases, these are obviously species with limited tolerance with respect to the environmental conditions in the country. Of particular interest are the rare species confined to the mountains. It can be supposed that they will be most vulnerable to the climate change. It can be expected that the ecological amplitude for these species will continue to shrink, and their abundance will decrease during climate warming. Rare species associated with caves will obviously be vulnerable to human disturbances in caves and their surroundings.

The beta coefficients for vulnerability (**Table 2**) show that in the lower and drier parts of the country, as well as near caves, there is a higher relative share of vulnerable species. As has already been mentioned, the anthropogenic negative impact on habitat characteristics important to bats is most pronounced in the lower parts of the country. The negative relation of LPC3 to species richness (**Table 2**) can be interpreted in this direction. The relative share of rare species is higher in sites covered with vegetation (lower mean scores of LPC3) and low in most places of anthropogenic origin (higher mean scores). Overall, the relative share of vulnerable species is greatest in the lower parts of the country, near caves, and in drier places with a relatively well-preserved natural vegetation.

[5] Peshev T, Peshev D, Popov V. Fauna bulgarica. Mammalia. Sofia: Editio Academica "Marin

Bats in Bulgaria: Patterns of Species Distribution, Richness, Rarity, and Vulnerability Derived…

http://dx.doi.org/10.5772/intechopen.73623

59

[6] Popov V. Terrestrial mammals of Bulgaria: Zoogeographical and ecological patterns of distribution. In: Fet V, Popov A, editors. Biogeography and Ecology of Bulgaria. Dordrecht:

[7] Johnson CN. Species extinction and the relationship between distribution and abun-

[8] Kerr JT, Kharouba HM, Currie DJ. The macroecological contribution to global change

[9] Myers N, Mittermeier RA, Mittermeier CG, da Fonseca AB, and Kent J. Biodiversity hotspots

[10] Boitani L, Maiorano L, Baisero D, Falcucci A, Visconti P, Rondinini C. What spatial data do we need to develop global mammal conservation strategies? Philosophical Transactions of Royal Society B Biological Sciences. 2011;**366**:2623-2632 http://dx.doi.org/10.1098/

[11] Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geo-

[12] Elith J, Graham CH, Anderson RP, Dudík M, Ferrier S, Guisan A, Hijmans RJ, Huettmann F, Leathwick JR, Lehmann A, Li J, Lohmann LG, Loiselle BA, Manion G, Moritz C, Nakamura M, Nakazawa, Overton JM, Peterson AT, Phillips SJ, Richardson K, Scachetti-Pereira R, Schapire RE, Sobeŕon J, Williams S, Wisz MS, Zimmermann NE. Novel methods improve prediction of species' distributions from occurrence data. Ecography. 2006;**29**:129-151

[13] Hernández PA, Graham CH, Master LL, Albert DL. The effect of sample size and species characteristics on performance of different species distribution modeling methods.

[14] Phillips S, Dudik M. Modeling of species distributions with Maxent: New extensions

[15] Peterson AT, Papeş M, Eaton M. Transferability and model evaluation in ecological niche

[16] Feilhauer H, He KS, Rocchini D. Modeling species distribution using niche-based proxies derived from composite bioclimatic variables and MODIS NDVI. Remote Sensing.

[17] Parviainen M, Zimmermann N, Heikkinen R, Luoto M. Using unclassified continuous remote sensing data to improve distribution models of red-listed plant species.

[18] Reports: Mammals. Available from: http://natura2000.moew.government.bg/Home/

modeling: A comparison of GARP and Maxent. Ecography. 2007;**30**:550-560

Drinov"; 2004. 632 p. (In Bulgarian)

dance. Nature. 1998;**394**:272-274

solutions. Science. 2007;**316**:1581-1584

for conservation priorities. Nature 2000;**403**:853-858

graphic distributions. Ecological Modelling. 2006;**190**:231-259

and a comprehensive evaluation. Ecography. 2008;**31**:161-175

Springer; 2007. pp. 9-37

rstb.2011.0117

Ecography. 2006;**29**:773-785

2012;**4**:2057-2075. DOI: 10.3390/rs4072057

Biodiversity and Conservation. 2013;**22**:1731-1754

Reports?reportType=Mammals [Accessed: 2017-12-03]

In this study, we quantitatively analyzed the Bulgaria-wide, high-resolution patterns in distribution of 29 bat species, species richness, rarity, and vulnerability. It was shown that species distribution models can effectively be used to reveal these patterns, covering areas that never have been sampled. The analyses revealed the individual role of important ecogeographic factors such as altitude, presence of caves, topographic wetness index, and various land cover features derived from satellite imagery. The results quantitatively confirmed the previously recognized types of distributional patterns, which were based on informal expert opinions. For the first time, high-resolution maps of species richness, rarity, and vulnerability were made.
