Revealing the role of past and current climate in shaping the distribution of two parapatric European bats, Myotis daubentonii and M. capaccinii

Abstract

In this contribution, we aim at determining the role of climate in shaping the current and historical (Last Glacial Maximum and Middle Holocene) distributions of two parapatric bat species (Myotis daubentonii and Myotis capaccinii) in Europe, using ensemble Habitat Suitability Modelling (HSM). Model projection to current climatic conditions predicted distributions largely coherent with the ones observed for the two species in the study area. Temperature- and precipitation-linked variables well accounted for the observed parapatry of the two target species. Moreover, areas of co-occurrence turned out to be those where the main ecological needs concerning the most important range-limiting variables are met for both species. Model projections to past scenarios allowed us to hypothesize the effects of climatic oscillations over the distribution of the two species during the Late Pleistocene and Holocene. Extreme range contractions during the Last Glacial Maximum and the subsequent expansions during Middle Holocene were predicted, consistently with general patterns of post-glacial recolonization routes of many temperate bat species in Europe. Our findings are largely coherent with recent phylogeographic studies investigating the two target species, thus corroborating the validity of potential distribution scenarios obtained from the models and, more importantly, confirming the great influence that glacial cycles had in shaping the biogeography of the European fauna.

Keywords:

• Myotis
• parapatry
• habitat suitability modelling
• past climate
• haplotype network

Introduction

Climate is often regarded as the main abiotic factor shaping species’ distributions both at broad (Root 1988; Buckley & Jetz 2007; Sexton et al. 2009; Wiens 2011; De Simone et al. 2020) and narrow (Iannella et al. 2018; Brunetti et al. 2019) spatial scales. Big climatic events, such as those that occurred during the Late Pleistocene and Holocene ages, have deeply influenced the distribution patterns of many animal and plant species, leading to large-scale range contractions, range expansions, and extinctions. Therefore, very often, current biogeographic patterns can be explained by analyzing species’ responses to those past climatic oscillations (Graham et al. 1996; Taberlet et al. 1998; Hewitt 2000; Stewart & Lister 2001; Prescott et al. 2012; Palkopoulou et al. 2013).

The effects of climatic oscillations, which took place during the Late Pleistocene and Holocene periods in Europe, are still detectable in many species’ current distribution patterns, and have been studied through fossil analysis, phylogeographic studies and palaeoecological reconstructions (Taberlet et al. 1998; Hewitt 1999; Sommer & Nadachowski 2006; Sommer & Zachos 2009; Alexandri et al. 2012). However, the distribution of a species in a certain time and space is determined not only by environmental factors but also by its interactions with other extant species (Case & Taper 2000; Peterson & Soberón 2012). In this context, parapatric species, which are characterized by adjacent ranges with limited overlapping areas, are often regarded as a model system for studying range limits and the factors shaping them. Many studies investigated the relative role of climate and biotic interactions in defining biogeographic boundaries of parapatric species (Huey 1979; Case & Taper 2000; Lawing et al. 2012; Werner et al. 2013; Cunningham et al. 2016; Iannella et al. 2017; Braz et al. 2019). When range limits are mainly shaped by abiotic factors, species tend to occupy all their environmentally suitable areas (“ecological parapatry” sensu Key (1981)), except for those geographically inaccessible (“migration” factor of the Biotic-Abiotic-Migration (BAM) diagram, Soberon and Peterson (2005)); usually, this type of parapatry is associated with relatively strong environmental gradients (Khimoun et al. 2013). Conversely, when range limits are mainly constrained by biotic factors (e.g. competition), environmentally suitable areas do not fully match the observed distributions and usually one species prevents the other from occupying common suitable areas (“competitive parapatry” sensu Haffer (1986)); in these cases, range margins can form at shallow environmental gradients (Case et al. 2005; Bridle & Vines 2007).

Habitat Suitability Modelling (HSM) techniques are among the most widely used tools for analyzing species’ distribution related to climatic conditions and, together with phylogeographic analyses based on genetic evidence, they allow a better comprehension of the processes determining species’ evolutionary history and distribution in time (Richards et al. 2007; Elith & Leathwick 2009; Araújo et al. 2019). Based on target species’ occurrence localities and environmental predictors, an HSM calibrates and infers the species’ response to these variables; then, it can project the calibrated scenario to new environmental conditions, which may vary in space and time. By allowing the identification of environmentally suitable areas for multiple species, HSM techniques have been proven to be very useful tools for studying parapatric range limits (Glor & Warren 2011; Khimoun et al. 2013; Iannella et al. 2017, 2018). Studies applying HSMs to investigate the factors shaping parapatric range boundaries have been conducted on many animal species, including mammals, birds, reptiles and amphibians (Cunningham et al. 2009; Gentilli et al. 2011; Glor & Warren 2011; Acevedo et al. 2012; Engler et al. 2013; Werner et al. 2013; McQuillan & Rice 2015; Iannella et al. 2017, 2018).

Despite their important ecological role (Ghanem & Voigt 2012), parapatry has rarely been studied in bats, and this is even more true if we consider HSM implementation in this context (Davis 1973; Baagøe 1984; Haffner & Stutz 1985; Rutishauser et al. 2012; Santos et al. 2014). Bats are particularly responsive to climatic and microclimatic changes in their environment as temperature- and precipitation-linked variables deeply affect every stage of their life cycle (Ransome 1971; Ibáñez 1997; Park et al. 2000; Kerth et al. 2001; Boyles 2007; Wojciechowski et al. 2007; Ben-Hamo et al. 2013).

In the present study, we assess the role of climate in shaping the current distributions of two parapatric bat species, Myotis daubentonii (Kuhl) and M. capaccinii (Bonaparte), using ensemble modelling (EM) techniques. These two vespertilionids show widely complementary distributions along the latitudinal gradient in Europe, with sympatric overlapping areas at their range limits.

For both target species, we build HSMs under the current climatic conditions using EM techniques to subsequently infer to past climatic conditions, namely the Last Glacial Maximum and the mid-Holocene scenarios, investigating whether the climatic oscillations that occurred in these periods might have influenced the two Myotis species, leading to shape their current distribution. Finally, we couple the HSMs’ outcomes with the current knowledge about genetic diversity aiming at understanding links between the current biogeographic and genetic assets for both species.

Materials and methods

Target species and study area

Two insectivorous bats, Myotis daubentonii (Kuhl, 1817) and M. capaccinii (Bonaparte, 1837), were considered in our analysis. They are among the so-called “trawling bats”, which hunt over waterbodies gleaning insects from the surface. These two species show wide complementary distributions in Europe, with overlapping areas at the edge of their ranges (Figure 1(a)) (Hutson et al. 2010; Kruskop et al. 2020).

Figure 1. Myotis daubentonii and M. capaccinii European distribution following the IUCN Red List data (Kruskop et al. 2020; Hutson et al. 2010) (a), their occurrence localities in the study area (b), the three-dimensional plot showing the overlap of their Hutchinsonian climatic niche and (c) binarized predictions for the current scenario, also considering the areas of sympatry (d)

M. daubentonii is widespread in central-northern Europe, whereas in the Mediterranean region, it displays a discontinuous distribution, mainly confined to mountainous areas. Since the environmental factors which might have affected this species in the past, conditioning the accessible area (sensu Barve et al. 2011), are explored here for the first time, only the Western Palearctic part of its distribution was considered. This species roosts in summer inside tree hollows, caves, buildings and other artificial structures, while its winter roosts comprise a wide range of underground habitats (Dietz & Kiefer 2016).

Conversely, M. capaccinii is a strictly troglophilous species, whose distribution is centered in the Mediterranean basin (Dietz & Kiefer 2016).

The study area is delimited by the coordinates 71.43 N-31.70 N, 10.00W-45.03E, covering each species’ known distribution in Europe and North Africa. This territory corresponds to the one considered by Rebelo et al. (2010), except for a slight extension southward (added to better encompass Myotis capaccinii’s distribution range).

Species data and climatic variables

Species’ occurrence localities were collected from the online databases GBIF and iNaturalist and integrated by conducting extensive bibliographic research for the whole study area, namely Western Europe and the Mediterranean basin.

All occurrences showing a coarse spatial resolution, equal or more than an area of 10 × 10 km, were discarded (Sillero & Barbosa 2020); further, to assess possible spatial autocorrelation among records, a Moran’s I test was performed in ArcMap 10.0 (ESRI 2010).

Bioclimatic variables of current, Mid-Holocene (MOL, ~ 6,000 years ago), and Last Glacial Maximum (LGM, ~ 22,000 years ago) periods were downloaded from the WorldClim.org website (Hijmans et al. 2005), version 1.4, at 2.5 arc-minutes resolution. Given their ecological relevance for the target species (Rebelo et al. 2010), we selected the following five bioclimatic variables: BIO1 (Annual Mean Temperature, °C*10), BIO2 (Mean Diurnal Range, °C*10), BIO6 (Minimum Temperature of Coldest Month, °C*10), BIO12 (Annual Precipitation, mm), and BIO17 (Precipitation of Driest Quarter, mm). To assess whether possible correlations among these predictors occur, a correlation matrix was built in ArcMap 10.0.

Two Global Climate Models (GCM) were considered for the Mid-Holocene and Last Glacial Maximum scenarios, the CCSM4 (Gent et al. 2011) and the MIROC-ESM (Watanabe et al. 2011), considering the uncertainties which could result from a single GCM approach (Marmion et al. 2009; Stralberg et al. 2015). Predictors’ values for current climatic conditions were extracted using the occurrence localities for both species through ArcMap 10.0 (ESRI 2010).

Niche overlap

To detect a possible climatic niche overlap between the two target species, we used the “hyperoverlap” package (Brown et al. 2020) in R (R Core Team 2016), developed on the basis of the n-dimensional Hutchinsonian niche concepts (Hutchinson 1957), and recently applied to biogeographic research (Iannella et al. 2021). The values of the aforementioned predictors, read on the occurrence localities of the two target species, were used as inputs for both the “hyperoverlap_detect” (kernel = polynomial; kernel.order = 3), and “hyperoverlap_lda” functions, to, respectively, quantify and display the possible niche overlap. The latter function builds a three-dimensional plot based on the linear discriminant analysis and the first two axes of the principal components’ residuals obtained by the former.

Model building and evaluation

To build the EMs, we used the ‘biomod2’ R package (Thuiller et al. 2016) in R (R Core Team 2016). The “ensemble” approach enables the use of different base modeling algorithms and averages the predictions of each resulting model based on the respective performance, obtaining more accurate and reliable outcomes (Thuiller 2003; Guisan et al. 2017; Cerasoli et al. 2019a; Iannella et al. 2019a, 2019b).

We generated 10 datasets with 1000 pseudo-absences each using the function ‘BIOMOD_FormatingData’ and “disk” strategy. This algorithm allows pseudo-absences to be randomly selected inside a predefined range of distance from occurrence points; in our case, a minimum distance of 100 km and a maximum distance of 1000 km were chosen (Barbet‐Massin et al. 2012; Cerasoli et al. 2019b; Thuiller et al. 2009), based on the spatial behavior of the target species (Barova & Streit 2018).

Three algorithms were selected to build the ensemble HSMs: Generalized Linear Models, Gradient Boosting Machine, and Random Forest, to encompass and take advantage of the different statistical approaches, from regression and classification to machine learning ones (Guisan et al. 2017). Every single model’s parametrization is reported in Supplementary file S1. Generalized Linear Models are among the most applied modelling techniques in HSMs (Thuiller 2003), while Gradient Boosting Machines (Friedman 2001) and Random Forests (Breiman 2001) are machine learning methods.

The predictive accuracy of the ensemble models was evaluated using both the area under the curve (AUC) of the receiver operating cCharacteristic (ROC) curve (Hanley & McNeil 1982) and True Skill Statistic (TSS) (Allouche et al. 2006). A 10-fold cross-validation was conducted, with 80% of occurrence records being used as training data and 20% as testing data. For both species, the percent contribution of each predictor was assessed, and the corresponding partial dependence plots were produced. The ensemble models (EMs) were then created from the single base models using the “BIOMOD_EnsembleModeling” function, selecting single models which performed with AUC > 0.80 and TSS > 0.70, thus obtaining a more “conservative” final EM (Iannella et al. 2018; Cerasoli et al. 2019a; De Simone et al. 2020).

Calibrated EMs were then projected to the LGM and MOL scenarios, by using both the CCSM4 and MIROC-ESM. For these past projections, the Multivariate Environmental Surface Similarity (MESS) (Elith et al. 2010) was calculated through the “mess” function of the “dismo” R package (Hijmans & Elith 2016), to obtain maps of the environmental dissimilarity (i.e. models’ extrapolation) of the past projections with respect to the calibration (current) scenario. These maps were then used to calculate the Multivariate Environmental Dissimilarity Index (MEDI) (Iannella et al. 2017), an algorithm that downweighs each GCM projection in the final averaged model based on the degree of extrapolation.

Habitat suitability maps were created using ArcMap 10.0 (ESRI 2010); a threshold value for the predicted suitability was calculated through the maximization of the TSS by the “ecospat.max.tss” function of the “ecospat” R package (Di Cola et al. 2017) to outline the two species’ binary distributions for the current scenario. These binarized maps were processed, along with the genetic information obtained from bibliographic sources (Viglino (2012) for M. capaccinii and Nardone (2015) for M. daubentonii), in GIS environment, which has shown to be a useful tool for several biodiversity issues (Iannella etal.2016, 2020; Rodríguez-Rodríguez et al. 2019; Rudke et al. 2020). To this aim, the predictions about haplotype network areas were obtained using the “Split binary SDM by input clade relationship” tool of the SDMtoolbox for ArcMap 10.0 (Brown et al. 2017).

Results

A total of 2558 occurrence localities for Myotis daubentonii and 423 occurrences for M. capaccinii were gathered (Figure 1(b)); these are, respectively, the 86.1% and 82.4% of the total occurrences obtained by applying a spatial resolution filter (Supplementary Material (Table S1)); no duplicate records were observed. Presence data showed no spatial autocorrelation (Moran’s I = −0.000813, p‐value = 0.395888, z‐score = −0.848988 for Myotis daubentonii and Moran’s I = −0.001402, p‐value = 0.699519, z‐score = 0.385970 for M. capaccinii); variables showed low-to-medium correlation values (max |r| observed = 0.74) (Supplementary Material (Table S2)).

The niche analysis (14% misclassified points, 51 Support Vectors) detected an overlap between the two species, showing a clear nestedness of the points distributed in the three-dimensional space (Figure 1(c)).

The models obtained showed good predictive accuracy, with AUC = 0.990 and TSS = 0.888 for M. daubentonii, and AUC = 0.977 and TSS = 0.832 for M. capaccinii. On the other hand, the “ecospat.max.tss” algorithm resulted in TSS_MAX = 0.454 for Myotis daubentonii and TSS_MAX = 0.467 for M. capaccinii, values used as thresholds to binarize (Figure 1(d)) the continuous models (Figure 2(a) and Figure 2b).

Figure 2. Predicted suitable areas for current climatic conditions for Myotis daubentonii (a) and M. capaccinii (b); distribution of the four shared and most contributing current bioclimatic variables for both species (c–f)

Maps resulting from the projection of the models to the current climatic conditions (Figure 2(a) and (Figure 2b)) report wide complementary areas, in which each species shows high levels of predicted suitability, as well as “shared” areas, with favorable bioclimatic conditions for both (Figure 2(a) and (Figure 2b)). These areas of sympatry mainly fall in the Balkans south of the Carpathian Mountains, in Italy, and along the southern coasts of France.

When comparing the binarized suitable areas (Figure 1(d)) with occurrence localities, the 55.1% (179) of the total 325 occurrences falling within the sympatric areas belong to M. capaccinii, while the 44.9% (146) belong to M. daubentonii. Also, this overlapping area represents 13.4% of the total (binarized) distribution of M. daubentonii and 46.1% of the one of M. capaccinii.

The relative contributions of the predictors used to build the models are reported in Table I for both target species, with some of them showing opposite and possibly discriminating trends in the response curves for the two Myotis species (Table I and Figure 3). In addition, these plots report the predicted suitability densities measured in the occurrence points falling within the sympatric areas: high predicted suitability in these areas is associated with intermediate values of the predictors (Figure 3).

Table I. Variables’ relative contribution inferred from the Ensemble Modelling process for the two target bat species, with BIO1 (Annual Mean Temperature), BIO2 (Mean Diurnal Range), BIO6 (Minimum Temperature of Coldest Month), BIO12 (Annual Precipitation), BIO17 (Precipitation of Driest Quarter)

Myotis daubentonii		Myotis capaccinii
Variable	Contribution (%)	Variable	Contribution (%)
BIO6	44.1	BIO1	38.1
BIO17	31.0	BIO12	23.2
BIO1	12.7	BIO6	19.0
BIO2	9.0	BIO17	13.6
BIO12	3.3	BIO2	6.1

Figure 3. Marginal response curves of the five predictors used for the modelling process and the density of predicted suitability read within the sympatric area in the occurrence localities, for Myotis daubentonii and M. capaccinii.

The habitat suitability map for the current scenario and past time-frames (MOL and LGM) are shown in Figures 2, 4 and 5, respectively; to better visualize and compare these trends with the most contributing predictors, we report their respective values for the study area in the corresponding side-panels (c–f) for each. It is observed that the most contributing predictors show marked variation especially at the LGM–MOL transition (Figures 4 and 5), with an apparent association with the contraction and subsequent expansion patterns of predicted suitability modelled for both species.

Figure 4. Predicted suitable areas for the Mid-Holocene climatic conditions for Myotis daubentonii (a) and M. capaccinii (b); distribution of the four shared and most contributing Mid-Holocene bioclimatic variables for both species (c–f)

Figure 5. Predicted suitable areas for the Last Glacial Maximum climatic conditions for Myotis daubentonii (a) and M. capaccinii (b); distribution of the four shared and most contributing Last Glacial Maximum bioclimatic variables for both species (c–f)

Both temperature- and precipitation-linked variables well account for the observed parapatry of the two target species, defining their range limits (Figures 2, 4 and 5). Meaningfully, the response curves for BIO1 (annual mean temperatures) display an opposite and “selective” pattern, with high suitability values relating to different temperature ranges for the two species. Likewise, the response curves for BIO2 (daily temperature range) reveal that areas of high predicted suitability for M. daubentonii tend to show a smaller diurnal temperature variation, which is typical of temperate-humid climates as defined by Köppen (Köppen 1900; Kottek et al. 2006), whereas areas of high predicted suitability for M. capaccinii show greater diurnal temperature variation, which is instead typical of the Mediterranean climate.

Concerning the precipitation-linked variables, M. capaccinii displays high values of predicted suitability in areas where annual precipitations (BIO12) reach at least 500 mm (Figure 3). Observing the BIO12 map for the current climatic conditions (Figure 2(e)), it is evident that mean annual precipitations are rather homogeneous in the whole study area, except for continental areas of Turkey and Spain and some coastal areas of southern Spain and northern Morocco. There, mean annual precipitations (and the corresponding suitability) display lower values; indeed, M. capaccinii has a more fragmented distribution in the western Mediterranean and is virtually absent in central Turkey, in line with what was reported from Dietz and Kiefer (2016). Moreover, when observing the predicted suitability and the precipitation of the driest quarter (BIO17), a pattern in the Mediterranean area can be observed for M. daubentonii: environmental suitability appears to be greatly limited where summers are very dry, a condition which is counterbalanced by other variables in mainland Europe (Figure 2).

The projections of the ensemble models to past climate conditions were made under the assumption of past niche conservatism (Wiens & Graham 2005), which has been already tested for bats (Rebelo et al. 2012; Razgour et al. 2013). These models allow us to hypothesize the effects of climatic oscillations on the distribution of the two Myotis target species during the Late Pleistocene and Holocene.

The Last Glacial Maximum projection reveals a strong decrease of the suitable areas for M. daubentonii, which resulted limited to the south-western portion of its current range, where the most range-limiting variables (i.e. BIO6 and BIO17) display possible appropriate values for the species (Figure 5(a, d and f)). This is due to enhanced evaporation over the North Atlantic, which in turn caused increased precipitations over southwestern Europe during the LGM (Ludwig et al. 2016). Conversely, very low temperatures and scarce precipitations in central-northern Europe, south of the Fennoscandian ice sheet, created adverse conditions to the presence of many temperate species in those areas (Ludwig et al. 2016). The predicted suitability for M. daubentonii was also found to be very low in the central and southern parts of the Iberian Peninsula, apparently due to the scarce precipitations of the driest quarter. Suitable areas are detectable also in Italy and the Balkans (Figure 5(a)), where temperatures and precipitations of the driest quarter display medium-to-high values (Figure 5(c, d and f)), thanks to a regional southerly atmospheric circulation which brought humidity from the Mediterranean Sea towards the Italian peninsula during the LGM (Giraudi 2011). Some exceptions are referred to the mountainous chains of the central Apennines and Dinaric Alps, where during the Last Glacial Maximum, relevant glaciers formed (Hughes & Woodward 2008).

The LGM projection for M. capaccinii reveals environmentally suitable areas limited to the southernmost coastal portions of its current range and to the Mediterranean islands, which are characterized by relatively high mean temperatures and medium values of annual precipitations (Figure 5(b, c and e)).

A certain degree of complementarity between the suitable areas detected for each species emerges in the LGM projections as well, coherently with variations observed for the most range-limiting variables. Furthermore, possible areas of sympatry emerge in southern Italy and the Balkans (Figure 5(a and b)).

The Middle Holocene projections reveal estimated distributions very similar to the current ones for both species (Figure 4(a and b)). Myotis capaccinii’s projections display high suitability areas slightly more extended, especially in the Iberian Peninsula and along the northern coasts of Africa, where monthly precipitations were more abundant than the current ones (Figure 4(e)).

The binarized models obtained for current climatic conditions (Figure 1(d)) were used to infer the haplotype network for both species (Figure 6). With regard to M. daubentonii, a trend of diversity that clearly decreases along the latitudinal gradient south-to-north is found, being especially high in the Italian and Iberian peninsulas (Figure 6(a)). Alongside this, the network for M. capaccinii (Figure 6(b)) reports the highest density of Voronoi polygons in the Balkans at the Turkish/Greek and Bulgarian borders. This density gradually decreases moving towards the west, with one vast polygon covering the binarized suitable areas in the Iberian Peninsula.

Figure 6. Haplotype network inferred from the binarized current predictions and the genetic information from Nardone (2015) for Myotis daubentonii (a) and Viglino (2012) for Myotis capaccinii (b)

Discussion

Models calibrated for both species showed good performances, with projections to current climate revealing estimated distributions that are largely coherent with the observed distributions of the two species in the study area. Our analysis of the EMs obtained for M. capaccinii and M. daubentonii clearly shows that their current distribution is strongly influenced by both temperature and precipitation-linked variables. This is in line with the responsivity of bats to climatic and microclimatic changes in their environment, as temperature- and precipitation-linked variables deeply affect duration and quality of their hibernation (from 4 to 15°C for M. capaccinii and from 0°C to 10°C for M. daubentonii) and daily torpor, awakening frequency, roost selection, and gestation length (Ransome 1971; Ibáñez 1997; Park et al. 2000; Kerth et al. 2001; Boyles 2007; Wojciechowski et al. 2007; Salari & Kotsakis 2011; Ben-Hamo et al. 2013). Overall, it appears that M. capaccinii’s distribution is related to the typical Mediterranean climate with hot-to-warm and dry summers, while M. daubentonii’s distribution is related to the temperate-humid climate characterized by warm-to-cool summers and weak seasonality in precipitations (Kottek et al. 2006). However, M. capaccinii’s range is primarily limited by mean annual temperature, and less by BIO17, being the species present also in Mediterranean areas where precipitations are more evenly distributed throughout the year and, therefore, summers less dry (as also resulted from the variables’ importance Table I). On the other hand, M. daubentonii is mainly limited by precipitations of the driest quarter, less by mean annual temperatures, being the species present also in Mediterranean areas with medium-to-high annual temperatures. Thus, common suitable areas turn out to be those where the main requirements relative to the most range-limiting variables are met for both species, high enough annual mean temperatures to ensure M. capaccinii’s presence and, at the same time, wet enough summers for M. daubentonii.

Despite being largely coherent with the observed distributions, the models display some discrepancies, particularly in certain regions of the Iberian Peninsula. In fact, M. daubentonii is marginally present in some areas of southern Spain, where our model predicts very low suitability values; conversely, M. capaccinii is virtually absent in Portugal and inland Andalucia, where our projections display areas of high suitability for this species. As a matter of fact, M. capaccinii’s distribution “anomalies” in the Iberian Peninsula have been already addressed by Garrido-García (2019), who hypothesized that the species’ absence in these environmentally suitable areas is related to changes in its historical distribution, between the Late Pleistocene and Holocene periods (see below).

As previously stated, parapatry between closely related species can be determined by abiotic factors, biotic interactions, or by a combination of both. In this respect, Biscardi et al. (2007) revealed a considerable niche overlap between our two target species due to similar foraging behavior and strategy, diet and, to a certain degree, ecological preferences. However, the question about whether the climatic niche overlap occurring between M. capaccinii and M. daubentonii results in competitive exclusion is still open; in fact, a partial overlap resulted from both HSMs and hypervolume-based analyses, as also observed in other cave-dwelling species in Western Europe (Mammola & Isaia 2017; Mammola et al. 2021). Compared to M. capaccinii, M. daubentonii is less selective with regards to roost choice, water quality and riparian vegetation’s conditions, making it a potentially stronger competitor. However, the high trophic niche overlap might not lead to competitive exclusion in sympatric areas because their main prey, Chironomids, can be considered an unlimited food resource (Krüger et al. 2012). Furthermore, their trophic niches may differ on a finer scale (e.g. distinct predation on adults and pupae, spatial segregation), allowing shared occurrence in the same area as documented by Krüger et al. (2012) for M. daubentonii and another trawling bat species, M. dasycneme Boie.

Given the relatively wide areas of sympatry observed for the two target species, and given the high degree of correspondence between the model-predicted distributions and the observed distributions, we can affirm that the species’ range limits in the study area are likely to be determined by climatic factors. This is also corroborated by the suitability densities measured in occurrence localities found in the sympatric area, graphically shown in Figure 3. It can be observed that, when responses to environmental predictors for the two target species intersect in a specific range of values, suitability reaches high values for both, as found from case studies on other taxa (e.g. Mammola & Isaia 2017; Iannella et al. 2017, 2018) and historical, macroecological trends (Hewitt 1999, 2000).

About the projections to the past climate, even though past niche conservatism was reported for some bat species (Rebelo et al. 2012; Razgour et al. 2013), some plasticity could likely have occurred, thus adding some hidden bias to our predictions.

Our results are coherent with the “refugium theory,” according to which many temperate species survived the glacial maximum in southern refugia and later re-colonized mid-latitude and Northern Europe from the Iberian Peninsula, Italy and the Balkans as confirmed for many other European bat species (Hulva et al. 2004; Juste et al. 2004; Rossiter et al. 2007; Rebelo et al. 2012; Razgour et al. 2013; Salicini et al. 2013). Considering the LGM and MOL projections, similar recolonization routes can be hypothesized also for Myotis daubentonii and M. capaccinii. Apart from southern refugia, our results allowed us to identify potential refugial areas in southwestern France, where LGM projections reveal highly suitable areas for M. daubentonii (Figure 5(a)). This hypothesis is also supported by the Late Pleistocene fossil remains found in the Dordogne region (Sevilla & Chaline 2011). As for M. capaccinii, we could identify potential refugial areas in the south-western portion of the Iberian Peninsula, northern Africa and along the southern coasts of Greece and Anatolia (Figure 5(b)).

Our results are coherent with the preliminary study carried out by Nardone (2015), which confirms that Mediterranean Peninsulas (Italy, Iberia, Balkan) acted as glacial refugia for M. daubentonii and suggests that European populations derived from the postglacial expansions of the Italian and Balkan lineages, while the Iberian lineage remained confined south of the Pyrenees. This is further confirmed by the haplotype network analysis, which reveals a higher haplotype diversity in the Mediterranean peninsulas (Figure 6(a)). Moreover, the density of the Voronoi polygons mirrors the predicted suitability levels of the LGM period in these regions: low to medium suitability levels are found in the Balkans, which harbor two haplotypes; high suitability levels are found in the Iberian and especially in the Italian peninsula, where the highest density of Voronoi polygons occurs (Figures 5(a) and 6(a)). Central and Northern Europe are almost entirely covered by the ‘h10’ Voronoi polygon which belongs to the Italian lineage, indicating a rapid post-glacial expansion of this northern-Italian haplotype. The ‘h2’ haplotype polygon, which covers the majority of the Balkan Peninsula and Eastern Europe, covers a wide area in Central Europe as well, suggesting its expansion to these regions; however, the presence of this haplotype in the United Kingdom should be further investigated, given that no samples were considered in the research of Nardone (2015). The Iberian Peninsula hosts a high polygon density, including also the ‘h10’ Italian haplotype and the ‘h2’ Balkan haplotype; however, when looking at the map, it is clear that the strictly Iberian haplotype polygons are found only in this region, indicating that no post-glacial expansion beyond the Pyrenees took place (Nardone 2015) (Figure 6(a)).

About M. capaccinii, molecular analyses carried out by Bilgin et al. (2008) allowed the identification of two glacial differentiation areas in south-eastern Europe, namely the Balkans and Anatolia, corresponding to two distinct haplogroups. The study carried out by Viglino (2012) adds to these haplogroups a third one grouping Sardinian, Spanish and South Italian haplotypes; according to the known M. capaccinii’s fossil records (Garrido-García 2019), this common western-Mediterranean haplogroup could be explained by the species’ disappearance in the Iberian Peninsula during the last glacial phase and subsequent recolonization from the Italian peninsula during the Holocene. A support to this last hypothesis results from the suitability predicted by our model on the Iberian Mediterranean coast, which resulted as low during the LGM and medium to high during the MOL period (Figures 4(b) and 5(b)). Furthermore, this would also explain M. capaccinii’s current absence from potentially suitable areas in the western part of the Iberian Peninsula (Garrido-García 2019).

M. capaccinii’s inferred glacial and post-glacial history is also corroborated as a likely scenario by the haplotype network analysis: the highest density of Voronoi polygons (corresponding to higher haplotype diversity) is found in the Balkans at the Turkish/Greek and Bulgarian borders, where a suture zone of at least two refugia was identified by Bilgin et al. (2008) (Figure 6(b)). Considering the LGM map (Figure 5(b)), it is likely that these refugia correspond to the southern coasts of Greece and Anatolia, where our models display two distinct areas of high predicted suitability. Fairly high haplotype diversity is found in Italy as well, with a large area in the central and southern parts of the country (including Sicily) harboring haplotype ‘h5’ and Sardinia and Corsica sharing the same haplotype polygon with the northern areas of Algeria and Tunisia. Viglino (2012) hypothesized that the presence of these haplotypes in Italy, both belonging to a unique haplogroup, is due to the movement of individuals from northern Africa to Sardinia and southern Italy, eased by the lowering of sea levels during glaciations. Although samples from northern Africa were not available in Viglino (2012) and, thus, further investigations are to be made, this hypothesis is likely to be true considering the high levels of predicted suitability displayed by our LGM projections in northern Tunisia and Algeria, which suggests the association of a strong refugium with the region. Such connectivity between Italy and the Ibero-Maghrebian region has been detected for several bat species (Bogdanowicz et al. 2015).

The Iberian peninsula is entirely covered by the ‘h5’ Voronoi polygon, which is present in central and southern Italy too. The very low haplotype diversity observed in this region would suggest a late post-glacial recolonization, in contrast with both our model projections and Viglino (2012), who identified a potential refugial area in southern Spain, but in agreement with Garrido-García (2019) according to which lack of fossil records dating back to the last glacial phase and the current absence of the species from potentially suitable areas in the western part of the peninsula suggests an Italian post-glacial recolonization (Figure 6(b)).

Conclusions

The influence of climate in shaping the distribution of our two target species emerges by the present study, also reinforced by the trends observed in the overlap area. The analysis of these shared areas was confirmed to be of primary importance when dealing with parapatric species, how evidenced also in other researches regarding several case studies on different taxa. The relevance of precipitation- and temperature-related variables was assessed in the co-occurrence areas, with analogous or opposite trends supporting the climate-driven distribution hypothesis, also confirmed by the genetic data. The climate-based HSMs’ current-to-past projections reported biogeographic assets which agree with the general glacial and post-glacial trends found for many species other than bats, thus giving new insights for the interpretation of the current distributions for this taxon. Complementary approaches confirmed our findings, giving interesting insights about the concurrent use and implementation of techniques (Peterson 2009; Pahad et al. 2020).

Further detailed phylogeographical studies are needed to comprehensively assess the two species’ evolutionary history in the study area; however, considering the acknowledged reliability of HSMs in this field, our results may provide a useful base for future investigations.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

Ministero dell’Istruzione, dell’Università e della Ricerca - A.I.M. Project PON R &I 2014-2020 No. 1870582.