Introduction

Climate is a significant driver of species distribution worldwide, with observed climate change driving shifts in geographic ranges. This is most visible in species that are specialists or have low adaptability to environmental changes1,2,3,4,5. Climate warming is more noticeable in the northern hemisphere6,7 as shorter and milder winters with decreasing number of days with snow cover and hot and dry summers8. Relatively high temperatures in winter have a negative impact on hibernating animals9,10, such as bats11. However, adverse conditions throughout the rest of the year can affect the availability of food and suitable shelter locations12.

Species distribution modelling techniques (SDM) are an effective method of predicting the impact of environmental change on species. SDMs use mathematical algorithms, along with species occurrence data and environmental variables, to predict ecological niches of species across both time and space13,14. Recently, the combination of SDM with nature conservation has emerged as a promising approach, offering valuable applications to understand the effects of environmental change on protected species15,16,17,18,19,20.

Among the various available modelling platforms, BIOMOD2 is widely used as an ensemble modelling framework to analyse species distributions13. By integrating multiple modelling algorithms, it helps mitigate methodological uncertainties, thereby increasing the robustness of predictions. SDMs employ a range of statistical and machine-learning methods to analyse species-environment relationships, including regression techniques (Generalized Linear Models—GLM21, Multivariate Adaptive Regression Splines—MARS22), machine-learning and complex algorithms (Artificial Neural Networks—ANN23, Generalized Additive Models—GAM24, Generalized Boosted Models—GBM25, Maximum Entropy—MaxEnt26, Random Forest—RF27, XGBoost—eXtreme Gradient Boosting Training28), classification methods (Classification Tree Analysis—CTA29, Flexible Discriminant Analysis—FDA30), and range envelope approaches (Surface Range Envelope—SRE31).

The environmental variables used in SDM are typically derived from climate datasets, with commonly used sources including CHELSA32,33, WorldClim34, and TerraClimate35. Such datasets differ primarily in spatial resolution, temporal coverage, and the types of environmental variables they provide, which influence their suitability for different modelling applications.

For analyses assessing the potential impacts of future climate change, General Circulation Models (GCMs) are essential36. These complex numerical models simulate the Earth’s climate system by integrating interactions between the atmosphere, ocean, land surface, and ice-covered regions. Developed by climate research institutions worldwide, GCMs form the basis for projecting future climatic conditions under different scenarios of greenhouse gas emission. Each GCM varies in its sensitivity to climate variables, which influences species distribution projections, and to reduce uncertainty, researchers often use ensembles of multiple GCMs rather than relying on a single model37,38,39,40.

The western barbastelle Barbastella barbastellus (Schreber 1774) is a psychrophilic species11,41,42, has a specialised feeding ecology, rendering it vulnerable to ecosystem change43,44,45,46,47. Its population is currently declining, and in 2023, it was classified as Vulnerable in Europe under the IUCN Red List of Threatened Species (criteria A2c)48. The species is also listed in Annexes II and IV of the European Union Habitats Directive. Consequently, it is crucial to protect not only existing habitats but also areas that may become suitable for the species in the future.

Previous work has predicted changes in B. barbastellus range estimates at the regional or country scale49,50,51,52. We extended this work by using species distribution models (SDM) to investigate the potential habitat suitability of B. barbastellus on the continental scale in Europe under current conditions and three future SSP scenarios: SSP126, SSP370, SSP585 (respectively low, intermediate and very high greenhouse gas emissions). Additionally, we accessed how climate change will affect potential habitat suitability in the areas where there is currently a network of protected areas Natura 2000 that are aimed at conserving this species.

Our research proposes three hypotheses. First, we hypothesise that the potential distribution of B. barbastellus in Europe will change in the future, with a noticeable shift toward northern regions. Second, we hypothesise that these changes will have negative consequences, resulting in significant reductions in currently suitable habitats. Third, we hypothesise that climate change will diminish habitat suitability within the Natura 2000 network of protected areas, which are currently designated to protect this species.

Materials and methods

Target species

The western barbastelle is an insectivorous bat, which belonging to the Vespertilionidae family. The characteristic black colour of its fur with silver or gold hair tips and the ears fused at the base (Fig. 1a,b) allow it to be clearly distinguished from other species42,53.

This species is strongly associated with forests, particularly old mixed or deciduous stands54,55,56. For roosts (in case of daytime and nursery roots) it prefers tree cavities, cracks, and places beneath peeling bark, favouring dead or partially dead trees (Fig. 1c,d); including deciduous species (oak, beech) and conifers (pine, spruce)57,58,59. Barbastella barbastellus hibernates in cold, poorly insulated natural and artificial underground spaces such as caves, shelters and wells11,60,61, as well as beneath tree bark62. It is among the first bat species to end hibernation, with its activity closely linked to the availability of early-flowering trees, such as willows, which attract moths46.

The western barbastelle specialises in hunting nocturnal Lepidoptera, which make up 90–99% of its diet43,45,63. Its unique echolocation, made through both the mouth and nose64, emits low-amplitude calls that allow it to effectively hunt moths with tympanate organs (e.g., owlet moths Noctuidae), which can detect and evade the echolocation of other bat species65,66,67. However, this echolocation strategy results in a short detection range, limiting its ability to hunt in open spaces. Instead, it relies on tree lines, hedges, and other linear landscape elements during migration to avoid exposed areas68.

Fig. 1
figure 1

Barbastella barbastellus with visible characteristic ears fused above the forehead (a) and in the top view (b), an example habitat—dead spruces due to the bark beetle outbreak in the Białowieża Forest (Poland) with protruding bark (c) and close up at roost under bark (d). The photos were taken by: Olga Łuczak (a, b) and Andrew Carr (c, d).

Full size image

Species occurrence data

We used a comprehensive distribution dataset for B. barbastellus, sourced from the Global Biodiversity Information Facility (GBIF) database69, as the basis for our analysis. Before integrating these data into our modelling process, a series of filtering steps were applied to ensure the quality and accuracy of the occurrences. Given the high-resolution nature of the climatic layers, we used (30 arc seconds, ~ 1 km), only occurrence locations with spatial accuracy exceeding this resolution were retained for further analysis. To address potential issues of model overfitting and to reduce environmental bias, we consolidated multiple species occurrences that fell within the same raster cell. These points were merged into a single representative occurrence, located in the centroid of the cell. Environmental filtering, which has been shown to outperform traditional geographical filtering approaches70,71, was implemented using an algorithm developed by Varela et al. (2014) available at https://github.com/SaraVarela/envSample. This process used principal component analysis (PCA) maps based on selected environmental variables (see Subsection Climatic layers) as filters. Species occurrence data points that fell within one unit of PCA1 and PCA2 were consolidated into a single randomly selected point. Additionally, we generated background data for the entire study area by converting raster pixels into point data using the ‘Raster to Point’ tool within QGIS version 3.34.572. As a result of these pre-processing steps, our final dataset comprised 1145 records documenting the occurrences of B. barbastellus across the study area (Fig. 2).

Fig. 2
figure 2

Spatial distribution of the final data set of Barbastella barbastellus occurrences used in the species distribution modelling. The map includes a scale bar in the bottom left corner, representing a distance of 500 km, and a north arrow in the upper right corner. The map was projected using the WGS 84 / Pseudo-Mercator (EPSG: 3857) coordinate system and generated using QGIS version 3.34.572.

Full size image

Climatic layers

We used the global dataset CHELSA V2.132,33, which provides bioclimatic variables with a spatial resolution of 30 arc seconds (~ 1 km), available online at www.envidat.ch. To ensure the independence of the variables and reduce the risk of multicollinearity, we applied the Pearson’s correlation coefficient to evaluate the relationships between the variables. Following a systematic exclusion process, any variables with a correlation coefficient of |r| ≥ 0.8073,74 were excluded. This filtering resulted in the selection of 11 bioclimatic variables, which were used in the species distribution modelling (Table 1).

Table 1 Bioclimatic variables selected for species distribution modelling of Barbastella barbastellus.
Full size table

To project future scenarios, we used downscaled Coupled Model Intercomparison Project Phase 6 (CMIP6) climatology datasets, sourced from the CHELSA V2.1 database32,33. Recognising the crucial role that General Circulation Models (GCM) play in climate modelling36, we adopted a strategy to reduce the uncertainties inherent in individual GCMs and strengthen the reliability of our projections. To address potential biases from selecting a single GCM, we applied an ensemble approach by integrating multiple models, including GFDL-ESM476, UKESM1-0-LL77, MPI-ESM1-2-HR78, IPSL-CM6A-LR79 and MRI-ESM2-080, and created averaged maps for each variable in specific Shared Socioeconomic Pathway (SSP) scenarios, providing a more robust analysis of future potential climate conditions. We focused on the final set of bioclimatic variables retained after correlation-based filtering (Table 1) and conducted projections that extended into the long-term future (ca. 2100). Within this period, we examined three distinct SSPs representing different levels of anthropogenic greenhouse gas emissions and their corresponding radiative forcing: SSP126 (low emissions), SSP370 (moderate emissions), and SSP585 (high emissions). These SSPs align with Representative Concentration Pathways (RCPs) 2.6, 7, and 8.5, respectively, offering a range of possible future climate trajectories.

Modelling procedure

To analyse both the present and future distribution of B. barbastellus, we applied an ensemble modelling approach using the BIOMOD2 package81. Since the modelling algorithms require both presence and absence data, we utilised the BIOMOD2::BIOMOD_FormatingData() function to generate 10,000 randomly selected background points, which were treated as pseudo-absences. To address potential bias arising from the random selection of pseudo-absences, this process was repeated three times82. Before selection, we excluded any locations with missing values for explanatory variables to ensure data integrity.

For the model calibration and evaluation, we employed the BIOMOD2::BIOMOD_Modeling() function. We implemented ten different algorithms: Artificial Neural Networks (ANN)23, Classification Tree Analysis (CTA)29, eXtreme Gradient Boosting Training (XGBoost)28, Generalized Additive Models (GAM)24, Generalized Boosted Models (GBM)25, Generalized Linear Model (GLM)21, Maximum Entropy (MaxEnt)26, Multivariate Adaptive Regression Splines (MARS)22, Random Forest (RF)27, and Surface Response Envelope (SRE)31. Each algorithm was configured with its default parameters. The occurrence data were divided into two subsets, 80% used for model calibration and 20% reserved for model testing. To minimise variability in model performance, we conducted five evaluation runs for each algorithm. Evaluation metrics, including the True Skill Statistic (TSS), which ranges from − 1 to + 1 and accounts for both omission and commission errors as well as the success of random predictions, and the Relative Operating Characteristic (ROC), which ranges from 0 to 1 and assesses model performance across all threshold values, were calculated using cross-validation. Higher TSS and ROC values indicated greater model performance82,83.

We used the BIOMOD2::BIOMOD_Projection() function to project calibrated models. In total, 150 models were generated, derived from three replicates of the pseudo-absence datasets, ten algorithms, and five evaluation rounds. To combine these outputs into a single ensemble prediction, we used the BIOMOD2::BIOMOD_EnsembleModeling() function. This ensemble model was constructed by selecting only the best performing individual models, which had a TSS score greater than 0.7, given the resistance of the TSS to the size of the validation set and the prevalence of the species82,84. The weighted mean approach was applied to combine these top models, thus excluding any underperforming ones and improving the overall performance of the ensemble model14,85.

The final step involved using the BIOMOD2::BIOMOD_EnsembleForecasting() function to generate ensemble projections for both present and future conditions, taking into account the three Shared Socioeconomic Pathway (SSP) scenarios. Additionally, we employed the BIOMOD2::get_variables_importance() function to quantify the contribution of each explanatory variable. The importance values were standardised to reflect their percentage contribution to the construction of the final ensemble model. All analyses were performed in R86 and the resulting maps were generated using QGIS version 3.34.572.

Spatial analyses

We used the BIOMOD2::BIOMOD_RangeSize() function to evaluate potential shifts in the geographical distribution of B. barbastellus under future climate scenarios. This function operates on the binary outputs of species distribution models, which are generated alongside ensemble projections using BIOMOD2::BIOMOD_EnsembleForecasting(). Binary projections are constructed automatically based on the evaluation metric value (in our case, TSS) with maximum specificity and sensitivity and then the corresponding suitability threshold is determined. Once the optimal threshold is identified (in our case, it was 530), the outputs of species distribution models (habitat suitability maps with pixel values ranging from 0 to 1000) are converted into binary maps, where pixels with habitat suitability values greater than 530 are designated as suitable (pixel value = 1), while those below this threshold are classified as unsuitable (pixel value = 0). To assess changes in species habitat suitability over time, BIOMOD2::BIOMOD_RangeSize() compares the number of currently suitable and unsuitable areas (pixels) with future projections. The function quantifies these changes into four categories: range loss (currently suitable habitats projected to become unsuitable; pixel value = − 2), stable presence (areas remaining suitable across both time periods; pixel value = − 1), stable absence (areas remaining unsuitable across both time periods; pixel value = 0) and range gain (currently unsuitable areas predicted to become suitable in the future; pixel value = 1). Additionally, BIOMOD2::BIOMOD_RangeSize() function performs a detailed pixel-by-pixel analysis, comparing the number of pixels classified as suitable or unsuitable in current projections with future projections. On the basis of this comparison, it calculates and summarises the percentage of species habitat gained and lost. For visualisation, we generated maps illustrating these projected range shifts using QGIS version 3.34.572.

Furthermore, we evaluated how future climate change may impact the potential habitat suitability of B. barbastellus within areas protected by the Natura 2000 network, where this species is currently conserved. First, we obtained detailed data from the European Environment Agency87, which provided a comprehensive list and shapefiles (representing the spatial distribution) of Natura 2000 sites where B. barbastellus occurs. We then cropped the outputs of species distribution models (habitat suitability maps with pixel values ranging from 0 to 1000) for both present and future conditions under the studied SSP scenarios to the boundaries of the Natura 2000 areas and extracted habitat suitability values for each pixel within a designated Natura 2000 site. Given the large dataset (n = 471,023), we performed an Anderson-Darling normality test using the nortest::ad.test function88 to evaluate the data distribution. Due to the non-normal distribution of the data, we applied paired Wilcoxon test to assess the significance of differences between present and future habitat suitability values in the specific Natura 2000 areas studied.

Results

Evaluation of algorithms performance and determination of climatic variable importance

Among the models used, GAM, GBM, CTA, and RF exhibited the highest performance, with both TSS and ROC values exceeding 0.8 (Fig. 3; Table 2). The remaining algorithms demonstrated weaker performance, with XGBoost and SRE showing the lowest predictive accuracy compared to the other models (Fig. 3; Table 2). The weighted mean approach used to construct ensemble projections ensures that models with higher accuracy are given greater weight, thereby enhancing the overall predictive capacity of B. barbastellus models.

Table 2 Evaluation of algorithm performance used in the species distribution modelling of Barbastella barbastellus, including the corresponding mean values and standard deviations (in parentheses) for the TSS and ROC evaluation metrics.
Full size table
Fig. 3
figure 3

Performance of the algorithm used in Barbastella barbastellus SDM, measured by ROC and TSS values. Dots represent the mean values, while the lines indicate the corresponding standard deviations for each evaluation metric. The figure was generated using R version 4.3.386.

Full size image

The distribution of B. barbastellus was predominantly shaped by two key climatic variables. Temperature seasonality (bio4) emerged as the most significant factor, accounting for 66% of the overall influence on species distribution (Fig. 4). Furthermore, the length of the growing season (gsl) was also an important factor, contributing 14% to the species’ distribution patterns (Fig. 4). Although the mean annual air temperature (bio1) had a relatively lower impact, it still contributed 5% to the distribution model (Fig. 4).

Fig. 4
figure 4

Percentage variable importance for the final ensemble model projection for Barbastella barbastellus. bio1 mean annual air temperature, bio2 mean diurnal air temperature range, bio3 isothermality, bio4 temperature seasonality, bio8 mean daily mean air temperatures of the wettest quarter, bio15 precipitation seasonality, bio18 mean monthly precipitation amount of the warmest quarter, bio19 mean monthly precipitation amount of the coldest quarter, gsl growing season length, gsp accumulated precipitation amount on growing season days, gst mean temperature of the growing season. The figure was generated using R version 4.3.386.

Full size image

Present and future potential distribution of Barbastella barbastellus with potential changes in its range

The current potential distribution of B. barbastellus is concentrated primarily in various regions of Europe, with the highest habitat suitability found mainly in the central and western parts of the continent (Fig. 5a). In southern Europe, the Iberian Peninsula is notable for its extensive high suitability, particularly in the coastal and northern mountain areas, spanning from the Northern Meseta through the Cantabrian Mountains to the Pyrenees, and extending south along the eastern coast (Fig. 5a). Similarly, the Apennine Peninsula, starting from the Alpine region, bypassing the Po Valley and extending along the Apennines down to Sicily, also shows significant habitat suitability (Fig. 5a). Furthermore, Corsica and Sardinia also exhibit high suitability (Fig. 5a). In Central Europe, the Armoricain Massif and the Massif Central in France, areas surrounding the main chain of the Alps, and the North European Plains display moderate to high suitability (Fig. 5a). The lowlands and coastal areas of Great Britain and Ireland also provide favourable conditions for the species (Fig. 5a). In the Baltic region, high suitability is predominantly observed in the southern parts of Sweden and the coastal regions of the Baltic states (Fig. 5a). On the contrary, eastern Europe presents a more varied pattern of suitability, with moderate suitability particularly noted in the Carpathian regions (Fig. 5a). The Balkan Peninsula also supports a favourable distribution, with substantial suitability in coastal and lowland areas around the Adriatic Sea, the Aegean coast, and the regions bordering the Black Sea, as well as the Caucasus Mountains (Fig. 5a).

Fig. 5
figure 5

Potential habitat suitability for Barbastella barbastellus in the studied area under current conditions (a) and three future SSP scenarios: SSP126 (b), SSP370 (c), and SSP585 (d). Darker colours indicate higher habitat suitability. Each map includes a scale bar in the bottom left corner, representing a distance of 500 km, and a north arrow in the upper right corner. The maps were projected using the WGS 84 / Pseudo-Mercator (EPSG: 3857) coordinate system and generated using R version 4.3.386 and QGIS version 3.34.572.

Full size image

Future projections indicate substantial changes in the distribution of suitable habitats for B. barbastellus throughout Europe under different climate scenarios (Figs. 5b–d and 6a–c). Under the SSP126 scenario (Figs. 5b and 6a), suitable habitats are predicted to expand northward, particularly in the Scandinavia and Baltic states. Central Europe is expected to remain largely suitable, while the Mediterranean regions may experience a slight decrease in habitat suitability (Figs. 5b and 6a). This scenario suggests a positive range change of + 13% (range gain + 28%, range loss − 15%), indicating a potential increase in suitable areas (Fig. 6a).

The SSP370 scenario, on the other hand, predicts a reduction in suitable habitats in most of the currently favourable regions, with particularly significant losses observed from the Mediterranean region through the Armorican Massif to the North European Plains, as well as in the Alpine region. However, it is projected that the most suitable habitats will persist in the lowland and coastal areas of Great Britain, Ireland and the Baltic region (Figs. 5c and 6b). There is a notable shift toward the northern parts of Europe, with increased habitat suitability in southern Scandinavia, northern Great Britain, and the Arctic Ocean region, especially along the coasts of the White Sea and the Barents Sea (Figs. 5c and 6b). This scenario results in a negative range change of -39% (range gain + 18%, range loss − 57%), reflecting a significant loss of suitable habitats (Fig. 5b).

Under the SSP585 scenario, most of the areas currently suitable are expected to become unsuitable, particularly in mainland Europe (Figs. 5d and 6c). Climate changes in this scenario are likely to cause there fragmentation of the remaining suitable habitats, with the most favourable areas for the species’ persistence becoming restricted mainly to mountainous regions, such as the Cantabrian Mountains and the Alps (Figs. 5d and 6c). The northern parts of the North European Plains, as well as southern Sweden and the British Isles, are projected to remain critically important for the species’ survival (Figs. 5d and 6c). This scenario also shows a shift towards the northern parts of Europe, with increased suitability in southern Scandinavia, northern Great Britain, and the Arctic Ocean region, particularly along the coasts of the White Sea and the Barents Sea (Figs. 5d and 6c). This scenario indicates a dramatic negative range change of − 56% (range gain + 17%, range loss − 73%), highlighting a severe reduction in suitable habitats across the species’ current range (Fig. 6c).

Fig. 6
figure 6

Changes in the potential range of Barbastella barbastellus under future conditions, as projected in the SSP126 (a), SSP370 (b), and SSP585 (c) scenarios, compared to present conditions. Each map includes a scale bar in the bottom left corner, representing a distance of 500 km, and a north arrow in the upper right corner. The maps were projected using the WGS 84 / Pseudo-Mercator (EPSG: 3857) coordinate system and generated using R version 4.3.386 and QGIS version 3.34.572.

Full size image

Changes in habitat suitability within Natura 2000 areas designated to protect Barbastella barbastellus

The analysis of potentially suitable habitats within Natura 2000 areas, where B. barbastellus is protected, aligns with general predictions for this species at the continental scale. The SSP126 scenario (low emission) indicates an increase in habitat suitability in these areas compared to current conditions (Fig. 7). In contrast, a significant decrease in habitat suitability is observed in the remaining future scenarios, with the decline becoming more pronounced as emissions increase (Fig. 7). The most substantial reduction occurs under the SSP585 scenario, representing a high-emission future, where habitat suitability is significantly lower than under current conditions (Fig. 7). This suggests that climate change, especially under higher-emission scenarios, will severely impact the species’ suitable habitats within protected areas, posing a potential challenge for conservation efforts.

Fig. 7
figure 7

Comparison of habitat suitability for Barbastella barbastellus under present conditions and three future SSP scenarios (SSP126, SSP370, SSP585) within areas conserved by the Natura 2000 network, where this species is currently protected. The boxplots illustrate the distribution of habitat suitability values for the present and future conditions. Each boxplot is paired with a density plot in the background. Habitat suitability is measured on the y axis, ranging from 0 to 1000. The black dot indicate the mean and the horizontal line the median. Significant differences between the present and future scenarios are marked with asterisks (****), indicating highly significant differences (p < 0.0001). The figure was generated using R version 4.3.386.

Full size image

Discussion

Among the climate-related variables shaping the potential distribution of B. barbastellus, those associated with temperature were the most influential, specifically, temperature seasonality, growing season length and mean annual air temperature. The first one emerged as the most critical factor to the overall projection of the species’ distribution model. This is significant because B. barbastellus is considered a psychrophilic species, particularly sensitive to temperature fluctuations11,41,42. The western barbastelle’s reliance on cold hibernation temperatures, typically between − 3.0 and 6.5 °C, with a preference for close to 0 °C, suggests that the temperature fluctuations are crucial for its survival41,42,89. Temperature seasonality might reflect the species’ preference for areas with stable microclimates, particularly for its underground hibernacula, where temperature depends on the average annual temperature in a given area12.

Recent studies show a division within the population into two wintering strategies. The more commonly observed strategy involves seeking colder, less insulated shelters11,90, which corresponds to the increase of the hibernating population in colder regions. This trend was reflected in our model and is associated with the change in the species’ range over time. The second wintering strategy, observed so far only in the ‘Nietoperek’ Natura 2000 site (Western Poland)—one of the largest bat hibernation sites in Europe—has been identified as a result of long-term research and is adopted by a minority of the population. It involves hibernating in warm (above 6.5°C), currently well-insulated hibernacula11. Hibernation at higher temperatures increases the metabolic rate, which accelerates the depletion of fat reserves, but also leads to more frequent awakenings. This strategy might be combined with foraging under favourable conditions. It should be noted, however, that our knowledge and understanding of barbastelle hibernation largely depend on the study site, the ability to identify and examine hibernacula, and the duration of observation in relation to the length of winter.

Barbastella barbastellus is one of the earliest to end hibernation in Central Europe91, emerging before many other bat species. As a result, the timing of its emergence coincides with the availability of early blooming plants, which are crucial for supporting insect populations. These insects serve as the initial food source for B. barbastellus after hibernation46. Therefore, the duration and timing of the growing season directly affect the availability of these essential food resources, making it a key factor for the survival and activity of the species in the early post-hibernation period. Deriving from the aforementioned ecological traits of the target species, we can identify the temperature seasonality as the most influential factor in our mode because it indicates the variability of temperature throughout the year, so the occurrence of seasons. The second most important to our model is the length of the growing season, which will influence the length of the period of food availability.

Projections suggest substantial range shifts for B. barbastellus in response to climate change. The northward expansion of suitable habitats aligns with range shifts observed in other bat species, such as Miniopterus schreibersii (Kuhl, 1817) and Pipistrellus nathusii (Keyserling & Blasius, 1839), which have also demonstrated shifts toward northern Europe in response to warming climates92,93. Similarly, increased habitat suitability in Scandinavia is consistent with general patterns of species migrating to cooler areas, where stable, cold hibernacula remain available94,95,96.

Predicted habitat fragmentation and significant loss of currently suitable habitats pose a major threat to B. barbastellus. This could result in severe population declines, particularly in southern Europe. This mirrors findings from studies on other thermally sensitive bat species, such as the temperate group represented by, among others, the western barbastelle, but most importantly the boreal group (Eptesicus nilssonii (Keyserling & Blasius, 1839), Myotis dasycneme (Boie, 1825), Nyctalus noctula (von Schreber, 1774), Vespertilio murinus (Linnaeus, 1758)), which are also expected to experience habitat losses due to climate-induced shifts97.

The western barbastelle is a priority species under the European Union Habitat Directive, and its protection through Natura 2000 sites underscores its conservation importance. Our findings indicate that as the range of the species shifts northwards due to climate change, conservation strategies must adapt to ensure the long-term survival of B. barbastellus. The predicted decrease in habitat suitability in a network of areas designed to protect this species further highlights this issue. In particular, regions projected to become suitable in northern Europe should be prioritised for future conservation efforts, especially those currently lacking in B. barbastellus populations. Furthermore, the maintenance of key landscape elements for this species, such as dead trees (also after bark beetle infestation59), mid-field tree stands58, rock with crevices on badlands with few trees98 or willow sites46 and early flowering trees and shrubs are especially important in areas where habitat fragmentation due to climate change could reduce the availability of essential roosting and foraging habitats.

Despite the robustness of the models used in this study, several limitations should be acknowledged. First, the models assume that the climatic niche of B. barbastellus will remain constant over time. This assumption neglects the potential for acclimatisation or adaptation to changing environmental conditions that may affect the future distribution of the species. Furthermore, although we were able to predict range changes based on climatic factors, the potential impact of changes in land use remains uncertain. Habitat fragmentation, driven by deforestation and other human activities such as agriculture and urbanization, along with the resulting loss of connectivity, further reduces the availability of suitable habitats for B. barbastellus. Moreover, the influence of new biotic interactions on the future distribution cannot be fully predicted. Further climate change with warmer winters and the choice of warmer shelters by part of the population may be determined to be more exposed to pathogens whose development has so far been limited by low temperatures in the hibernaculum, such as white nose syndrome (WNS). For example, in ‘Nietoperek’, the average temperature in the main corridor in the 1980s was around 8 °C99, but at present it has risen to 10 °C11. Low temperatures (under 5 °C) and high temperatures (above 19.0 °C) have been observed to limit the growth of the fungus Pseudogymnoascus destructans ((Blehert & Gargas) Minnis & D.L. Lindner, 2013) causing WNS100,101,102. It reaches the highest growth rate at 12.5–15.8 °C on media plates100; however, studies on bats in the wild101 and in the laboratory102 indicate a growth peak already at a temperature of 5–6 °C. Barbastelles can be exposed to completely new pathogens for this species, the spread of which will be contributed by thermophilic bat species expanding their range92,95,97,103. Additionally, western barbastelle could face competition for shelter access from an increasing number of representatives of opportunistic species in the future. A similar situation is currently observed in the daily barbastelle roosts occupied by Pipistrellus pygmaeus (Leach, 1825) in Wolin National Park104. Finally, our model does not account for molecular variation between populations, which may have an effect on their adaptability105.