Assessing future shifts in habitat suitability and connectivity to old-growth forests to support the conservation of the endangered giant noctule

Study area and spatial data

We selected the western part of the Palearctic region (longitude range: ~−10° ÷ 60° E) as the study area. This choice was determined by combining information about Nyctalus lasiopterus current range, as assessed within the IUCN Red List of Threatened Species (Alcaldé, Juste & Paunović, 2016), with occurrence data from published literature and online resources.

In detail, the species’ range spans from Portugal, Spain, France, Italy, Greece, Cyprus, and Turkey (with few scattered historical records also in northern Morocco and eastern Libya), to various countries in the Balkans (Slovenia, Croatia, Albania), central and eastern Europe (Slovakia, Poland, Belarus, Hungary, Bulgaria, Romania, Ukraine, European Russia), also extending to the Caucasus (Georgia, Azerbaijan) and western Kazakhstan (Alcaldé, Juste & Paunović, 2016).

The species’ occurrence localities were gathered through literature search and the Global Biodiversity Information Facility (GBIF, 2022). As for the literature data, we excluded old records (before the 1970s) and those lacking precise geographic information about the occurrence locality (e.g., articles generally reporting the presence of the species within a Region/Province or an entire Protected Area). Moreover, we excluded GBIF records being duplicates or reporting geographic coordinates with a spatial resolution coarser than 2.5 arc-min (i.e., ~5 km). These filters permitted to select N. lasyopterus occurrence records with a spatial and temporal resolution matching that of the predictors used to fit the ENMs (Sillero & Barbosa, 2020).

To assess the climate-occurrence relationships of N. lasiopterus, we downloaded a set of 19 raster layers representing bioclimatic variables from the Worldclim 2.1 (Fick & Hijmans, 2017) archive (https://www.worldclim.org/data/index.html, accessed on 9 March 2022), at 2.5 arc-min resolution for the ‘current’ timeframe (i.e., 1970–2000 average climatic conditions) as well as for three future time horizons (i.e., 2030, 2050, and 2070). For each future horizon, we downloaded raster data representing predicted climatic conditions under two Shared Socioeconomic Pathways (SSPs). Specifically, we selected the SSPs 3.70 and 5.85: the latter predicts an increase in Earth’s radiative forcing of 8.5 W m⁻² by 2100, due to lack of international effort in limiting greenhouse gas (GHG) emissions, and is frequently named ‘business as usual’ scenario (Riahi et al., 2017); differently, SSP 3.70 is usually labelled as ‘middle of the road’ because it predicts an increase in radiative forcing up to 7.0 W m⁻² in the next decades, in a context of “regional rivalry” with limited global cooperation to mitigate anthropogenic global warming (Riahi et al., 2017).

Predictions from the climate-based ENMs were then refined in a post-modelling phase (see below) by including topographic and habitat-related predictors. As for topography, we downloaded a Digital Elevation Model (DEM, 25 m resolution) from the European project Copernicus repository (https://land.copernicus.eu/imagery-in-situ/eu-dem, accessed on 9 March 2022). The habitat-related variables were instead derived from two distinct datasets: (i) vector spatial data about European rivers from Grill et al. (2019), from which we selected the major riverine areas (i.e., RIV_ORD > 4 based on the Strahler stream order, indicating the level of branching of watercourses in a top-down fashion) as these latter are used by N. lasiopterus to roost, forage, and move across vast landscapes (Popa-Lisseanu, Bontadina & Ibáñez, 2009); (ii) occurrence localities of European old-growth forests from Sabatini et al. (2018). Occurrence localities of old-growth forests were converted from point to polygon features through a three-step procedure: first, we derived the radius (in meters) of the considered forest from the corresponding area, reported in hectares in Sabatini et al. (2018), through the formula $\sqrt{\left(oldgrowthforestarea\ast 10000\right)/\pi }$; then, we used this radius to calculate a buffer around the occurrence point of the forest; finally, we refined the obtained polygon by discarding the areas which were encompassed by the buffer but did not fall into forests’ features of the European Nature Information System (EUNIS) spatial dataset (100 m resolution), which we downloaded from the European Environment Agency website (https://www.eea.europa.eu/data-and-maps/data/, accessed on 9 March 2022). Details about each of the considered old-growth forests, including the extent of the corresponding refined polygonal feature, is reported in Table S1.

Ecological niche modelling

The ecological niche modelling step was performed in R version 4.1.1 (R Core Team, 2021). Specifically, we took advantage of the ‘gbm’ R package ver. 2.1.8 (Greenwell, Boehmke & Cunningham, 2019) to fit the ENMs for N. lasiopterus based on the gradient boosting model (GBM) algorithm, known also as boosted regression trees (Elith & Graham, 2009). The choice of this algorithm derived from both the possibility of easily tuning it in the R environment and to the fact that it emerged as one of the best performing ENM algorithms with presence-pseudo-absence data, once properly tuned (Elith & Graham, 2009; Hao et al., 2020).

To feed the GBM algorithm, we selected a subset of the 19 Worldclim bioclimatic variables by integrating the results of a Variance Inflation Factor analysis (VIF), performed through the ‘usdm’ R package (Naimi, 2015), with relevant published information about the species’ ecology (Cerasoli et al., 2021). We chose VIF ≥ 10 as exclusion criterium for the single variables as it is deemed a suitable threshold to deal with multicollinearity in Ecological Niche Modelling (Guisan, Thuiller & Zimmermann, 2017). Then, we generated 5,000 pseudo-absences through the ‘disk’ strategy of the ‘BIOMOD_FormatingData’ function of the ‘biomod2’ R package (Thuiller, Georges & Engler, 2016), setting 130 and 250 km as the minimum and maximum radius, respectively. The rationale of this choice was threefold. First, selecting pseudo-absences within a geographically constrained buffer permits to avoid them falling either too close to the presence points, which could hamper a proper model calibration, or too far from them, which could artificially inflate model performance estimates (VanDerWal et al., 2009). Secondly, GBM was shown to perform well when pseudo-absences are located one or two degrees away (i.e., roughly 120-to-240 km within the latitudinal range of our study area) from presence points (Barbet-Massin et al., 2012). Lastly, including ecological knowledge about the dispersal capability of the target species usually decreases the risk of pseudo-absences being instead false absences (Phillips et al., 2009): as the maximum daily foraging transit recorded for N. lasiopterus is about 130 km (Naďo et al., 2019), selecting pseudo-absences outside this minimum distance from presence points reduces the risk of labelling as absences localities potentially used by the giant noctule populations comprised in our occurrence dataset.

Then, we weighted presences and pseudo-absences so that the sum of the weights of the former equals the one of the latter. Indeed, assigning a same overall weight to presences and pseudo-absences usually increases ENMs’ predictive performance when the generated pseudo-absences are far more numerous than the available presences (Cerasoli et al., 2017; Gouvêa et al., 2020; Thiault et al., 2020).

To detect the best parametrization for the GBM algorithm, we built three different matrices containing several combinations of ‘gbm’ parameters and the respective set of values (for brevity, here we list the ones for the first matrix only: shrinkage = 0.01, 0.1, 0.3; interaction.depth = 1, 3, 5; n.minobsinnode = 5, 10, 15; bag.fraction = 0.65, 0.8, 1). Then, we ran as many GBM models as the combinations, increasing the n.trees value from 10,000 to 20,000 but keeping the train.fraction = 0.8 and the cv.folds = 10 fixed. We finally chose the set of parameters resulting in the lowest root mean square error (RMSE) (Friedman, 2001; Greenwell, Boehmke & Cunningham, 2019; Cervellini et al., 2021). To investigate possible spatial autocorrelation in predictions from the optimized GBM model, which could have resulted for instance from sampling bias affecting N. lasiopterus occurrence data (Roberts et al., 2017), we fitted an additional set of 10 GBM models, using the tuned set of parameters. In each of the 10 iterations, the model was fitted on the occurrence data coupled with a set of 5,000 pseudo-absences randomly selected within the above-mentioned 130-to-250 km buffer. Then, we computed residuals on the calibration data for each of these models. Subsequently, we draw a correlogram, through the ‘ncf’ R package version 1.2-9 (Bjornstad, 2020), showing the variation of spatial autocorrelation, as represented through the Moran’s index (I), at increasing average inter-point distances. Finally, we computed the mean and median pairwise Euclidean distances across N. lasiopterus occurrence data to check whether they fall within the inter-point distance range not affected by spatial autocorrelation in model predictions, as emerging from the correlogram.

Successively, we assessed the discrimination power of the optimized GBM model through the Boyce index (Boyce et al., 2002), which is particularly suited for ENMs built on presence and pseudo-absence data (Hirzel et al., 2006; Leroy et al., 2018). Moreover, we measured the relative contribution of the selected variables through the randomization algorithm implemented in the ‘summary.gbm’ function of the ‘gbm’ R package.

Then, we projected the optimized GBM model across the entire study area for both current climatic conditions and various future scenarios represented by the combinations of year (2030, 2050, and 2070) and SSP (SSP3.70 and SSP5.85). We buffered possible differences in future predictions, which may be observed when comparing ENMs’ projections obtained using distinct General Circulation Models (Stralberg et al., 2015), by considering three of them, namely the BCC-CSM2-MR (Wu et al., 2019), the IPSL-CM6A-LR (Boucher et al., 2020) and the MIROC6 (Tatebe et al., 2019). Moreover, we first assessed the prediction uncertainty due to extrapolation (i.e., the divergence between environmental conditions at calibration points and those across the projection surface) through the Multivariate Environmental Similarity Surface (MESS) (Elith & Leathwick, 2009) implemented in the ‘dismo’ R package (Hijmans & Elith, 2016). Then, we used the resulting MESS maps to implement the Multivariate Environmental Dissimilarity Index (MEDI). This index weighs ENMs’ projections under different GCMs based on the corresponding MESS, finally returning a combined weighted projection (Iannella, Cerasoli & Biondi, 2017). We iterated this process of ENMs’ refinement for each year ^× SSP combination. ENMs’ predicted suitability, ranging from 0 to 1 (low-to-high suitability), was then reclassified on a 1-to-10 scale for post-modelling purposes.

Post-modelling analyses

To sharpen the ENMs-derived predicted suitability for N. lasiopterus, we applied the “couple-and-weigh” framework from Iannella et al. (2021a) by integrating spatial information about elevation and rivers’ network. All the post-modelling processes subsequently described were performed in ArcGIS Pro 2.9 (ESRI Inc, 2022).

First, we extracted elevation values at occurrence localities from the DEM so as to draw an elevation preference curve, normalizing the “raw” occurrence frequencies along the elevation gradient (100-m-wide intervals) to a 1-to-10 scale. Then, we took advantage of the information about “Connectivity Status Index (CSI)” and “Urbanization (URB)” present in the rivers’ dataset of Grill et al. (2019) to calculate a connectivity-functionality index (obtained as $CSI-URB$), in turn converted to a 1-to-10 scale. In this manner, spatial data about elevation and rivers’ network had the same scale of values (i.e., 1-to-10) as the reclassified predicted suitability from the climate-based ENM.

Finally, to obtain current and future weighted suitability maps, we combined the current and future (MEDI-corrected) predictions from the climate-based ENMs with the reclassified DEM and rivers layers through the ‘weighted overlay’. In this step, the cell resolution of the elevation and rivers layers was automatically upscaled to that of the ENM predictions (i.e., 2.5 arc-min) by the ‘weighted overlay’ tool itself.

Landscape connectivity

Starting from the obtained current and future weighted suitability maps, we calculated connectivity across the study area through the Circuitscape v.5 (Anantharaman et al., 2020) package in Julia programming language (Bezanson et al., 2017). Circuitscape, which was shown to perform well in ecological connectivity studies (Dickson et al., 2019), applies coupled random walk and circuit theory-based algorithms to model ecological connectivity across landscapes using a suitability (or friction) map as a base, and target/destination points (or areas) to be connected as nodes (McRae et al., 2008; McRae, Shah & Edelman, 2016). Two main node categories were used, with subsequent sub-algorithms applied, as we assessed both connectivity among N. lasiopterus occurrence localities and connectivity from these latter towards European old-growth forests. For the former aim, we applied the ‘pairwise’ approach so that, for each possible pair of occurrence localities, each term is, iteratively, the source or the destination node (McRae, Shah & Mohapatra, 2013). For the latter aim, we chose the ‘Advanced’ mode to select N. lasiopterus current populations as sources and old-growth forests as destination areas (McRae, Shah & Mohapatra, 2013). Then, we calculated, for both the pairwise and the directional (i.e., occurrences-to-forests) case, the Standardized Connectivity Change Index (SCCI) between current and future scenarios, which ranges from −1 (connectivity loss) to +1 (gain), with 0 representing the connectivity stability (Iannella et al., 2021a, 2021b).

Considering the daily distances that the giant noctule may cover (Naďo et al., 2019), the connectivity estimates were refined by applying a 130 km buffer around old-growth forests to highlight those falling within this distance from N. lasiopterus occurrence localities (i.e., the forests being easily reachable by the species).

Bioclimatic niche overlap and gap analysis

To assess possible climatic links between N. lasiopterus populations and old-growth forests, we used the ‘hyperoverlap’ R package (Brown, Holland & Jordan, 2020). Specifically, we segregated occurrence localities of N. lasiopterus and old-growth forests in three “longitude groups” (Longitude ranges: −15°–0°, Western; 0°–15°, Central; 15°–50° Eastern). Then, we first assessed overlap in bioclimatic space among the three groups of N. lasiopterus occurrences (42 Support Vectors, SVM kernel = linear), by adapting the ‘hyperoverlap_set’ function (Brown, Holland & Jordan, 2020) to perform the analysis with more than two groups (modified function reported in Supplemental Material). Successively, we assessed niche overlap between each longitude group of N. lasiopterus occurrences and each category of old-growth forests (e.g., ‘Near-virgin Forest’, ‘Long Untouched Forest’), divided by longitude group as well (104 Support Vectors, SVM kernel = polynomial, 2^nd degree).

Also, for each longitude group we measured: (i) the distance between the N. lasiopterus occurrence localities and the old-growth forests; (ii) the values of the most contributing bioclimatic variables (as resulting from ENMs) within the occurrences of both N. lasiopterus and old-growth forests.

We finally evaluated the protection status of N. lasiopterus by performing a gap analysis (Scott et al., 1993; Jennings, 2000; Hermoso et al., 2022). Specifically, we assessed the number of occurrence localities falling within protected areas, as well as their distance from the easily reachable old-growth forests (i.e., those located within the 130 km buffer). For this purpose, we downloaded spatial data about the Natura 2000 sites and Nationally Designated Protected Areas (https://ec.europa.eu/environment/nature/natura2000/data/index_en.htm and https://www.protectedplanet.net/en/thematic-areas/wdpa?tab=WDPA, respectively; accessed on 9 March 2022).

Occurrence data and climate-based ecological niche models

At the end of the filtering phase, we selected 260 occurrence localities for Nyctalus lasiopterus (Fig. 1, Table S2).

Average inter-point distance within the filtered occurrence dataset equaled 1,562 km, while the median inter-point distance was 1,349 km. The correlogram drawn on residuals from the 10 additional runs of GBM fitting indicated that spatial autocorrelation was practically absent (i.e., Moran’s I ~ 0) when inter-point mean distance ranged between 1,250 and 2,000 km (Fig. S1). Thus, we may be quite confident that no spatial bias possibly hidden in the occurrence dataset noticeably affected the ENM fitting process and the resulting predictions.

Based on the VIF scores (Table S3), the climatic variables selected for model fitting were BIO1 (annual mean temperature), BIO2 (mean diurnal temperature range), BIO6 (minimum temperature of the coldest month), BIO12 (annual precipitation), and BIO17 (precipitation of the driest quarter). The lowest RMSE was recorded for the GBM model fitted with 15,380 trees, interaction.depth = 3, shrinkage = 0.001 and n.minobsinnode = 15, which also attained Boyce index = 0.867. The three most contributing variables within this model were: BIO6 (38.2%), which showed a sigmoid-shaped curve of predicted suitability peaking at ~−1 °C; BIO1 (18.3%), with suitability linearly decreasing after ~0 °C and BIO2 (17.5%), showing a peak in the 6–12 °C interval (marginal response curves are reported in Figs. S2A–S2C). Also, BIO6 – BIO2 represented the pair of variables attaining the highest interaction score (0.385), with highest predicted suitability when BIO6 ≥ 5 °C and BIO2 ranges between 8 °C and 12 °C (Fig. S1D).

Weighted model

After the implementation of the “couple-and-weigh” step, in which the reclassified map from the optimized climate-based ENM was merged with the reclassified information about elevation and major rivers, the resulting weighted model predicted highly suitable areas (i.e., weighted suitability >0.8) within an extent spanning −9° ÷ 35° in longitude and 36° ÷ 47° in latitude (Fig. 1), and scored a Boyce index = 0.903.

Under current conditions, areas in the southern and northern Iberian Peninsula, southwestern France, peninsular Italy, coastal areas of the Balkans, and coastal Anatolia host the highest weighted suitability (Fig. 1). On the contrary, only small suitable patches are predicted in central and eastern Europe, corresponding to occurrence clusters.

A general decrease in weighted suitability on the gathered N. lasiopterus occurrence localities emerged from current to future scenarios, with suitability losses being more prominent as the time horizon of model projections moves forward (Figs. 2A and 2B). In 2030, the percent suitability loss compared to current conditions is essentially equal under the SSP 3.70 (−8.05%) and SSP 5.85 (−7.61%) scenarios. In 2050 and 2070, predicted losses are even larger, bringing the total reduction to −17.90% and −17.31%, for the SSP 3.70 and SSP 5.85 scenarios, respectively.

The western areas, in particular those in the range of longitude −8° ÷ 30° and latitude 36° ÷ 48°, will experience the broadest suitability losses, with stronger evidence under the SSP 5.85 scenario (Figs. 2C and 2D). Our models, on the other hand, predict a significant increase in suitability in the north-eastern areas, particularly in the range of longitude 26° ÷ 45° and latitude 44° ÷ 60°, especially under the SSP 5.85 scenario, most likely as a result of possible temperature increases in the northernmost latitudes (Figs. 2C and 2D).

Landscape connectivity among Nyctalus lasiopterus populations

The connectivity analysis among N. lasiopterus occurrence localities, based on the weighted suitability map obtained under current conditions, resulted in a dense, vast connection between northeastern Spain and southwestern France, overriding the Pyrenees (Fig. 3A). Moreover, this connectivity core is linked to other smaller corridors systems in southern France (Haut-Languedoc chain and coasts of southwestern Mediterranean France) and the Azahar Coast. On the other hand, smaller intra-connections emerged in Andalusia, Julian pre-Alps, and Bosphorus, even though there are few occurrence localities in the surroundings of the two latter. Also, the pre- and sub-Alpine Italian areas showed medium connectivity values, potentially linking the western systems to the eastern ones. On the contrary, the central European (Slovakian) occurrence cluster resulted in lower intra- and inter-connections.

When analysing future changes in connectivity through the SCCI, a loss can be observed in the southern part of the Iberian Peninsula–southwestern France corridor system, as well as in the corridors emerging in western Andalusia, Bosphorus and Slovakia (Fig. 3B). On the other hand, concurrent stability is predicted for the Haut-Languedoc chain and coasts of southwestern Mediterranean France. In contrast, connectivity increases in the northern part of the Julian pre-Alps system (Fig. 3B). Also, significant connectivity gains are observed for northern Portugal, southern Spain and northern France; in a similar way, central and eastern Europe gain connections across broad areas. All these mentioned trends are more pronounced as the time horizon for model projection is shifted farther in the future (i.e., connectivity changes are greater in 2070 than in 2030); also, the SSP5.85-based predictions generally resulted in higher connectivity changes than those based on SSP3.70, especially in areas showing currently high connectivity.

Connectivity towards old-growth forests, niche overlap, and protection status

About the connectivity inferred from the actual N. lasiopterus occurrence localities towards old-growth forests, we found that major corridors exist where old-growth forests are scarce. In fact, Sabatini et al. (2018) reported a high number, or even clusters, of old-growth forests in central Italy, the Balkans, and central-western Europe, while the major corridor systems emerging from our analyses are located in the northern Iberian Peninsula, southern France and Italian Alps (Fig. 4A). In addition, another small corridor system is observed in Andalusia, which is more linked to the farther old-growth forests in the north than to the nearer ones, located eastward.

The SCCIs calculated on connectivity modelled between N. lasiopterus occurrences and old-growth forests located within the 130-km buffer show overall corridors’ stability in the future except for the Pyrenees, where a slight increase in connectivity is predicted for all the scenarios (Fig. S3A), and for Andalusia, where instead a slight loss is predicted. Similarly, the areas outside the 130-km buffer show stable connectivity, except for the gains reported in the eastern Alps, northern France, and Denmark (Fig. S3B).

The three bioclimatic variables contributing the most to the climate-based ENM (i.e., BIO6, BIO1 and BIO2) show decreasing values along the west-east longitude gradient, when sampled in N. lasiopterus presence localities. A similar trend emerges for the same values measured in old-growth forests’ occurrences, altough these latter show higher variability (Fig. 4B). A mixed trend is instead observed for the distances between N. lasiopterus occurrences and old-growth forests, with the eastern ones being closer and the central and western ones occurring at more considerable distances (Fig. 4B).

In the bioclimatic space modelled through the ‘hyperoverlap’ R package, we found a separation (i.e., non-overlap) between the western and eastern N. lasiopterus populations, with the central ones overlapping with both (Fig. 5A). When examining bioclimatic niche overlap between the distinct longitude groups of N. lasiopterus and old-growth forests, eastern N. lasiopterus populations overlap with all the forests’ categories except the western ones, whereas the western N. lasiopterus populations overlap only with western forests (Fig. 5B). Central N. lasiopterus populations overlap with ‘Old-Growth Forest’ and ‘Long Untouched Forest’ categories from both the eastern and central longitude groups. Moreover, N. lasiopterus eastern populations are the only ones overlapping with the “near-virgin forest” category (Fig. 5B).

Finally, we found 59 N. lasiopterus occurrences (out of 260) falling into protected areas, of which 15 in Natura 2000 sites, eight in Nationally Designated Protected Areas, and 36 covered by both (i.e., where these two kinds of protected areas overlap). The median distance between these N. lasiopterus occurrence localities and old-growth forests ranges from 20 km to nearly 120 km, depending on the considered longitude group and category of protected areas (Fig. 6).