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

Abstract

Bats are affected by a variety of anthropogenic pressures, and effective conservation measures require a complex approach by not only covering the roosts themselves but also the surrounding habitats and migration corridors. The development of concrete localized conservation measures requires detailed quantitative data to assess habitat status regarding the most crucial factors for concrete species. This chapter aims, by modeling using a Maxent based on many georeferenced locations and the state of ecologically relevant ecogeographic variables, to reveal the spatial trends in the habitat suitability of 29 bat species; to obtain meaningful biogeographical species groups; and to provide a countrywide quantitative assessment of bat richness, rarity, and vulnerability. The modeling results showed that altitude, karstic areas (presence of caves), topographic wetness index, and presence of deciduous forests were the most influential factors. In this respect, three well-defined groups were delineated. The species’ richest areas were mostly located in semimountain karstic areas with a well-developed broadleaved forests, and the lowest in xerophilous, bare habitats, especially those of anthropogenic origin. Regarding rarity, more rare species were associated with caves and mountains. Vulnerability (in terms of IUCN criteria) was positively affected by the presence of caves showing the importance of protecting karstic areas.

1. Introduction

Bulgaria has a uniquely high diversity of bats. Of the 35 species that are found in continental Europe, 33 species are found in Bulgaria. This is largely due to the transitional geographic location of the country, the diverse habitats, the significant elevation gradient from the sea level to the altitude above 2900 m, the preserved wildlife in many parts of the country, especially in the mountainous and semi-mountainous regions, and the presence of over 5900 caves. All the bat species are strictly protected by the National Biodiversity Act (Annex 3). Twelve species are listed in Annex 2 of the Habitat Directive. Despite the legal protection, many species have conservation problems. The main threats identified both for bats as a whole and for individual species, fall into two general categories such as anthropogenic influence on the roosts and habitat loss and degradation [1].

The first group of threats mainly concerns the cave-dwelling species. Many of them are vulnerable to human impacts because they are often more visible and found in higher numbers. Anthropogenic disturbance and vandalism, excessive caving visits, destruction, placing doors and bars that prevent or hinder the partial or full access to them, agricultural and animal breeding activities in caves, and underground water catchments are among the major potential threats. Measures to mitigate the impact of these threats are directed at identifying important underground roosts and their legal protection. Out of the 5900 caves in Bulgaria, there are about 125 caves and cave’s complexes declared as Natural Monuments. Among them, 52 caves are known as Important Bat Underground Habitats of national and 13 of international importance [2]. Within the borders of 17 Protected Sites, there are at least 120 caves, many of them with importance to bats. Additionally, 817 caves fall within the borders of National and Nature Parks and 173 caves are part of Strict or Managed Nature Reserves. The most important bat underground sites according to Habitats Directive (92/43/EEC) are designed as Natura 2000 network sites. Some of the other Natura 2000 sites, particularly those covering large areas, also contain many bat caves [1].

Regardless of the formal legal protection of important bat caves and their immediate surroundings, it still affects only a small portion of the bat’s diversity. In fact, the fate of the protected areas and the biological diversity they contain is influenced to a great degree by actions within the surrounding landscape. Effective conservation measures require a complex approach, covering the habitats around caves, including bats’ foraging habitats and migration corridors.

Conservation of the habitat is also important for a large group of bats that are not directly related to caves and whose ecology is little known. Many human activities are potentially negative in this respect—destruction of natural vegetation especially forest, including clearing of the few remaining natural lowland forests for agricultural purposes and of older forests in the higher elevations for timber; widespread alternations of mid-elevation forests due to clearing, fires, heavy pressure from livestock grazing, and artificial planting and forestation (especially the replacement of broad-leaved forests with conifer plantations). Destruction of some habitat elements has a particularly negative impact on bats. Felling old tress and trees with hollows restrict the possibilities of finding appropriate roosts, especially for the nursery colonies; destruction of the natural open water areas (lakes, marshes, river arms); destruction of the hunting habitats and the flight corridors. Negative for bats are many activities associated with agricultural land uses: plowing of meadows, including formerly uncultivated lands; overgrazing; expansion of monocultures and input-intensive agriculture, especially the intensified use of fertilizers and pesticides; poorly planned construction and development projects, including wind turbine construction, tourist resorts and facilities, highways and other transportation projects, mines, and quarries, as well as urban expansion in general [1, 3].

The development of concrete measures for integrated bats conservation requires better insight into the environmental requirements of the species. It is necessary to identify locally specific measures. Given that there are differences in the ecological characteristics of the species in different regions, specific quantitative data are needed. The detection of distribution patterns along environmental gradients is an important task in conservation ecology. By knowing species-environment relationships, species and species assemblages can be used in understanding the conservation needs of poorly known species with a narrow niche breadth. Although such analyses and generalizations aimed at identifying groups of species with similar ecological requirements already exist [1, 4, 5, 6], such classifications were made by eye and were based entirely on expert judgment.

With the advent of increased interest in numerical classification, clustering of multivariate species data became very popular in such studies. To be effective, this approach needs to base on comprehensive quantitative data on the distribution of individual species. Such data, however, are not always available, especially for rare species and those with a hidden lifestyle such as most bat species. In recent years, habitat models relating habitat characteristics, in the form of digital coverage of ecogeographic variables, and species occurrences or abundances are increasingly used for estimating habitat suitability and forecasting species distribution. Moreover, this approach, based on niche theory, has proven useful in understanding the rules governing species assembly at various spatial scales. The search for causes determining patterns in species distributions in natural and disturbed landscapes is of primary importance in ecology, and establishing relationships between species distributions and environmental characteristics is a widely used approach. Modeling also plays an increasingly important role in conservation [7, 8], particularly for understanding impacts of global change on biological diversity, identifying gaps in protected area networks, and for planning and reserve design [9]. Furthermore, the model approach provides the opportunity to obtain high resolution maps that are particularly important for terrestrial conservation planning, where cell sizes of 1–100 km² are commonly required, depending on the organism and local habitat heterogeneity [10].

Recently, the author [1] modeled bat species listed in Annex 2 of the Habitat Directive across the country using a 0.63 km2 resolution. The study was based on location data with precise geographic coordinates available to date, mainly those published in the summary article of Benda et al. [4], using three modeling methods such as ecological niche factor analysis, generalized linear model, and discriminant analysis.

In recent years, more effective modeling methods based on presence-only data have become increasingly popular. Among these methods, Maximum Entropy, a recently developed modeling 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 autocorrelation.

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, compared with climatic/topographical variables-only models [16].

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 improved by including this type of data [17].

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 existing knowledge gaps.

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 bat conservation.

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 Directive.

2. Material and methods

2.1. Species data

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].

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.

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.

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.

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).

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].

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.

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 well-developed 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.

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.

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 km² or less [34, 35, 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

$$\mathrm{W}=exp{\left(-\frac{{\mathrm{Q}}_{\mathrm{i}}-{\mathrm{Q}}_{\min }}{r_i\times {\mathrm{Q}}_{\max }-{\mathrm{Q}}_{\min }}\times 0.97+1.05\right)}^2$$

where Q_i is the occurrence of species i, Q_min and Q_max 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

$$\mathit{\text{OIRR}}=\frac{\frac{\sum w_i}{S}-w_{\min }}{w_{\max }-w_{\min }}$$

where w_i is the weight of the i-th species in the assemblage, S is the assemblage species richness, w_min and w_max are the minimum and maximum weight, respectively.

Abundance-based IRR (AIRR) was calculated as

$$\mathit{\text{AIRR}}=\frac{\frac{\sum a_i{w}_i}{N}-w_{\min }}{w_{\max }-w_{\min }}$$

where a_i and w_i 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 w_min and w_max 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{{\displaystyle \sum _{i=1}^S}{v}_{\mathit{ri}}}{S_r}$$

where v_n is vulnerability rank of species i and S_r 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 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.

Species	N-locs	Niche breath	Training AUC	Test AUC	Alt	Cave	TWI	LPC2	LPC1	LPC3
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
Rhinolophus ferrumequinum	464	0.94	0.88	0.86	16.73	47.76	8.22	2.67	22.69	1.92
Rhinolophus hipposideros	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 ecogeographic variables to the Maxent model.

Contrary to expectations, there was no correlation between the niche breadths and the AUC-test 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

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.

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.

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.

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.

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).

Figure 5.

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

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.

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.

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.

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.

Figure 6.

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

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.

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 R^2	0.82	0.59	0.46

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.

4. Discussion

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.

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.

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.

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.