An iconic messenger of climate change? Predicting the range dynamics of the European Bee-eater (Merops apiaster)

When environmental conditions change, species usually face three options: adaptation, range shifts, or extinction. In the wake of climate change, it is generally believed that range shifts are the norm in mobile species such as birds, resulting in poleward range shifts. The European Bee-eater is a predominantly Mediterranean species which has expanded its range to higher latitudes over the last decades. Germany in particular has seen a surge in breeding pairs and foundation of new colonies. However, while many experts suggest climate warming as the main driver behind this range expansion, an explicit quantification remains open. Here, we use an ensemble modelling approach to study the recent climatic niche suitability of the European Bee-eater across Europe with a special focus on Germany and project its predicted Palaearctic breeding distribution onto the year 2050 using two global circulation models and two representative concentration pathways. Models were able to predict the current European range of the species with some underestimated areas in Central and Eastern Europe, depending on the selected model. We found a strong relationship between climatic suitable areas and estimated population sizes across European countries that is reflected in most algorithms. In particular, the German population size is in line with climate suitability in the country suggesting a strong climate–population relationship and a high degree of niche filling. Most future predictions point to an ongoing northward expansion of the species while areas in Southern Europe and the Maghreb area remain largely suitable. The strong climate–population relationship makes the European Bee-eater an appropriate indicator species for climate change. Yet the high variability of modelling algorithms also call for caution of using these techniques without careful inspection.


Introduction
Changes of environmental conditions make species usually react in three non-exclusive ways: evolutionary adaptation, range shifts or extinction. For species with potentially high dispersal abilities, such as birds, range shifts are believed to be the norm in the wake of climate change. As a consequence, breeding ranges of species in higher latitudes may shift polewards (Thomas and Lennon 1999;Parmesan and Yohe 2003;Hickling et al. 2006;Thomas 2010) or upwards along elevational gradients (Chen et al. 2011;Freeman et al. 2018). Both effects are tantamount with the assumption that species are able to track their climatic niches over time. Niche tracking can be observed over different time scales, ranging from intra-annual variation of seasonal niches up to niche conservatism across phylogenetic timescales (Peterson 2011;Engler et al. 2017). Nevertheless, while a growing number of studies support the assumption of niche tracking in response to climate change (e.g. Araújo et al. 2005;Tingley et al. 2009; La Sorte and Jetz 2012; Monahan and Tingley 2012), the question whether niche tracking is the norm is still unresolved . This is because even highly mobile species might be restricted in their dispersal capacities and track their niches only after a certain time lag (Devictor et al. 2008).
The European Bee-eater Merops apiaster is a turanianmediterranean faunal element which mainly breeds in the southwestern Palaearctic region in mediterranean, steppe and desert zones (Voous 1960;Fry and Boesman 2020). As long-distance migrants European Bee-eaters spend the northern winter in the Afrotropical region, mainly in western and southern Africa (e.g. Fry 2001). While the species also breeds in Namibia and South Africa (Fry 1988), we focus on the European breeding distribution of the species where it has shown a remarkable ongoing range expansion over the last decades at its northern range edge in Central and Western European countries such as Germany, Switzerland, France, Poland or Slovakia (e.g. Glutz von Blotzheim and Bauer 1980;Krištín and Kaňuch 2005;Gerber et al. 2011;Gedeon et al. 2014;Essel et al. 2016). Due to its excellent data coverage and its representative status as recently populated area, Germany will be in particular focus of our analysis. European Bee-eaters often breed in colonies and are generally inhabitants of open landscapes including agricultural land and river valleys (Voous 1960;Cramp 1985). The birds burrow nests in bare ground, often in banks, scarps and embankments, but also choose secondary habitats like quarries for breeding (Cramp 1985;Bauer et al. 2005). Most of these habitat types are locally exposed to warmer microclimates (e.g. Heneberg and Šimeček 2004).
While species' ranges are often related to climatic factors on a large scale (Pearson and Dawson 2003;Pigot et al. 2010), this could be particularly evident in the European Bee-eater. European bird species with a southern range centroid are often regarded as climate change indicators under global warming (Gregory et al. 2009). Moreover, trends of common bird species are related to their climate niche characteristics ) offering a role for precipitation and temperature as driver of population change in Mediterranean birds (Herrando et al. 2019). Indeed, current and past range shifts in the European Bee-eater have been related to climate change (Kinzelbach et al. 1997;Reif et al. 2010), and the species has been called a 'flagship species for a northward range extension of a southerly distributed bird species in Europe' (Fiedler 2016). However, an explicit quantification of the relationship between the current range expansion of the European Bee-eater and climatic factors issurprisingly-lacking. Even more so, other factors than climate might affect southern species which extend their ranges into Central Europe. For example, the range extension of the Melodious Warbler Hippolais polyglotta might at least partially be explained by biotic interactions with its sister species, the Icterine Warbler Hippolais icterina (Engler et al. 2013). Other southerly distributed species like the Great White Egret Ardea alba showed an increase in abundance in Central Europe, most likely due to a combined effect of habitat change, the cessation of persecution, conservation measures and climate change (Ławicki 2014). Disentangling the interplay of different factors for the European Bee-eater hence requires a detailed analysis of climate-occurrence relationships.
We apply a set of species distribution models (SDMs) within an ensemble approach to: (1) study the current climatic niche suitability of European Bee-eaters across their Western Palearctic breeding range and (2) project its breeding distribution onto presumed climatic conditions of the year 2050 using two global circulation models (GCM) and two representative concentration pathways (RCP) of greenhouse gases. Over the last two decades, species distribution models have become a major tool for forecasting ranges of birds and other taxa (Guisan and Thuiller 2005;Elith et al. 2006;Elith and Leathwick 2009;Engler et al. 2017), although their application for range shifting species remains challenging (Elith et al. 2010). A special case of SDMs is ensemble models in which possible uncertainties of single SDMs are counterbalanced by combining multiple algorithms (Thuiller 2003;Araújo and New 2007). Moreover, SDMs were successfully applied to infer population dynamics such as species abundance allowing to quantify possible functional relationships of species with their environment (VanDerWal et al 2009;Thuiller et al. 2014;see Engler et al. 2017 for an exhaustive overview in birds). Herein, we present and compare results of eight different SDM algorithms used to quantify current climate-population relationships and predict the future range of the European Bee-eater. Specifically, we analyse whether climate models cover the current breeding range of the species and whether areas of niche unfilling at the northern range edge are prevalent. We assume that current population sizes in Europe are closely related to climatically suitable areas. Our objective is to test if the population size of the European Bee-eater in Germany follows this assumed relationship with climatic suitability. If the prediction of the correlation is correct, the observations will be close to the fit (Fig. 1, case B). Alternatively, higher or lower suitable areas than expected by the current population size would point to factors different from climate as a main driver of population dynamics, such as niche expansion (Fig. 1, case A) or lack of habitat or food resources (Fig. 1, case C). By applying future predictions, we want to know whether the northward shift of the species' range continues and whether this is accompanied by a range contraction at the southern range of its distribution.

Species records
We used two sources of occurrence data to achieve a broad coverage of the palaearctic breeding range of the European Bee-eater Merops apiaster, the Global Biodiversity Information Facility (GBIF) and a database on distribution and population size of German breeding colonies. Specifically, we collected (1) georeferenced species data from GBIF (https:// doi. org/ 10. 15468/ dl. 4tjk09, last access 04.09.2019). We included "preserved specimens" and "human observations" from the core breeding season (i.e. June and July) to exclude most transients as well as vagrants out of their regular breeding area. To reduce sampling bias in presence records (Fourcade et al. 2014a), we filtered records using a grid of 10° by randomly picking one record per grid cell. Moreover, we excluded any data outside the time frame 1979-2013 so that environmental and species data are in congruence (see below). While we focus on the European breeding distribution of our species, we also include occurrences from the eastern part of the range to improve niche estimations (Barbet-Massin et al. 2010). We limit records to a rectangle from 18° western longitude to 100° eastern longitude and 20°-66° northern latitude. Note, that we only infer predictions for currently inhabited ecoregions further restricting our study area (see below). (2) Since the GBIF dataset showed some geographical gaps at the northern and western distribution in Central Europe, we took additional occurrences of breeding colonies from the extensive database of the "Special Interest Group European Bee-eater" of the German Ornithologists' Society (H.-V. and A. Bastian, unpublished data). To the best of our knowledge, these data represent a nearly complete coverage of the breeding distribution of European Bee-eaters in Western, Southern, and Central Europe. In sum, we used 1382 unique occurrences for further analyses (Fig. 2).

Fig. 1
We generally expect a strong positive correlation between suitable area (y-axis) and population size (x-axis) as the basis for the range expansion of European Bee-eaters with climate change (case B). If a region or country has strongly deviating residuals (cases A and C), European Bee-eaters follow other factors rather than climate 1 3

Environmental data
We used high-resolution climate data available through the CHELSA V1.2 database (www. chelsa-clima te. org; Karger et al. (2017a, b)). We restricted the climate data to the months May to September to cover the breeding season of the species in the Western Palearctic. Based on monthly data of this five-month period, we created six predictor variables which are largely analogous to the familiar so-called bioclimatic variables (Hijmans et al. 2005): minimum temperature, maximum temperature, average temperature, precipitation sum, precipitation of the wettest month, and precipitation of the driest month. Present climate covers the time period 1979-2013. For future predictions, we applied four climate change scenarios for the year 2050 based on two different widely used climate global circulation models (GCMs: CCSM4 and MIROC5) and two representative concentration pathways (RCP 4.5 and 8.5). While RCP 4.5 represents an intermediate assumption, RCP 8.5 can be described as a very high emission scenario (IPCC 2014).

Modelling
We applied an ensemble modelling approach using eight different SDM algorithms as implemented in the Bio-mod2 R package (based on Biomod,  including generalized linear models (GLM), generalized boosted models (GBM), random forests (RF), classification tree analysis (CTA), multivariate adaptive regression splines (MARS), artificial neural networks (ANN), flexible discriminant analysis (FDA), and the maximum entropy approach according to Phillips et al. (2006); i.e. MAX-ENT.Phillips; see  for details). As model parameterization needs to be within an ecologically meaningful context (Guisan et al. 2014), we generated 10,000 random background records in ecoregions covered by > 95% of European Bee-eater presence records during the breeding season (Fig. 2). We excluded ecoregions with fewer records, because the species breeds in river valleys in deserts or other specific microhabitats which are not well covered by climate data and are hence likely to bias predictions. We extracted ecoregion information from the WWF Terrestrial Ecoregion map (Olson et al. 2001, available through https:// www. world wildl ife. org/ publi catio ns/ terre strial-ecore gions-of-the-world). In particular, we used the 'randomPoints' function in the dismo R package (Hijmans et al. 2016) that allows for a latitudinal correction since the environmental layers are not in a planar coordinate system resulting in a projection bias towards higher latitudes (Hijmans et al. 2016). We ran each algorithm 10 times, randomly keeping 30% of presence records for model testing in each iteration. We Fig. 2 Overview of the study area with the background (dark grey) covering ecoregions incorporating 95% of species records (red) during the breeding period. To test the role of climate-driven north-ward range expansion, we focus on Germany (dark shade) where an exhaustive monitoring of breeding Bee-eaters is in place left all other model settings at default mode apart from rescaling predicted occurrence probabilities from 0-1 to a 0-1000 integer scale for memory saving. For model quality assessment, we used the Area Under the Curve (AUC) of the receiver operating curve (ROC, Fielding and Bell 1997) and true skill statistic (TSS, Allouche et al. 2006) classification as evaluation metrics. We ensembled models for each algorithm separately applying a TSS threshold of 0.5 and using committee averaging. This allowed us to quantify algorithm-specific differences while applying the same quality criteria. To restrict for possible extrapolation errors in the predictions, we clipped areas outside the ecoregions chosen as background (see Guisan et al. 2014 and above) and used MESS maps (Elith et al. 2010) to clip areas of non-analogous climate in the dismo R package. We finally converted predictions into binary presenceabsence maps, using the maximum sum of sensitivity and specificity thresholds among all algorithms as performed in Biomod2. We ultimately stacked the batch of maps for each scenario into consensus maps.

Country-level analyses
To elaborate the general relationship between the climatically suitable area and population size, we first compiled information on country-level estimates of breeding populations of European Bee-Eaters (BirdLife International 2015). In particular, we selected country-level information with 'good' or 'medium' quality ranks while discarding country-level information ranked as 'poor'. This way, we gathered information from 20 countries. For these countries, we calculated the number of suitable cells for each SDM algorithm (i.e. cells with suitability score > threshold) based on the current climate. Given that values in both predictors stretch across several orders of magnitude, we log(x + 1) transformed each before including them in linear models to assess the strength of relationship for each SDM algorithm. Using the predict() function in R, we estimated confidence intervals to check if the position of the German population in that relationship deviates strongly in either one or the other direction (Fig. 1).

Model performance
According to a summary of statistical parameters, all model algorithms performed well with AUC values above 0.8 and TSS values above 0.5 (Table 1). Minimum precipitation had a consistently high variable importance across all algorithms. Likewise, mean temperature was ranked high for six out of eight algorithms where it was also ranked highest among all predictors. In the two algorithms (i.e. GBM and MAXENT) that considered the mean temperature less important, minimum temperature became critical (Table 1).

Current potential distribution
The European distribution of the European Bee-eater is well covered by the stacked summary map, and population Table 1 Summary statistics of model performance (mean ± sd) for 8 algorithms across 10 iterations and their respective ensembles.
Shown are also the threshold values under Maximum Sensitivity and Specificity (MaxSS) as well as the normalized variable importance for each predictor and algorithm. For better visualization and comparability, importance values are highlighted by bars (blue = positive, red = negative, size = strength). strongholds of the species on the Iberian Peninsula, in Italy and Romania show high suitability values (Fig. 3). This pattern is obvious for all model algorithms, although there are some differences in suitability for Central and Northern Europe, mainly for the MAXENT and GBM algorithms, which mostly show unsuitable conditions north of the Alps (Fig. 3). In addition, outside of our focus area, large parts of the realized distribution in Asia and easternmost Europe were not recognized as suitable.
Focussing on Germany, all models were able to predict suitable areas along the upper Rhine valley. However, the breeding population in Eastern Germany was not predicted by all models (Fig. S1). Furthermore, patterns of overprediction are obvious in Northern Europe, often in coastal areas around the Baltic Sea and the North Sea.

Climate suitability and population size
Linear regression models showed significant positive relationships between the climatic suitable area and population size across European countries ( Table 2). The average coefficient of determination was considerable (R 2 = 0.411) and ranged between R 2 = 0.287 in MARS and R 2 = 0.653 in GBM (Table 2). For most algorithms, the German population size of the European Bee-eater falls within confidence limits of this relationship with climatic suitability (Table 2; Fig. 2). Only two algorithms, ANN and RF, overpredicted the amount of suitable area in Germany, and GBM slightly underpredicted the area (Table 2). A few other countries show strong deviations from these relationships, namely: Belarus (for CTA, FDA, GBM, GLM, MARS, MAXENT), Latvia (for CTA, GBM, MAXENT, RF), and Lithuania (for CTA, GBM and MAXENT), where the respective algorithms failed to predict suitable distribution and hence underestimated population sizes in that country (Fig. 4).

Future potential distribution
The RCP4.5 models predict a northward shift of the distribution of the European Bee-eater for both circulation models (Fig. 5a, b). Differences between MIROC and CCSM are small, with MIROC pointing to high suitability areas in Central and Eastern Europe for nearly all model algorithms. As with the current distribution, model algorithms also differ slightly in potential future distribution (Fig. S2). Following the RCP8.5 scenario, areas in the south of the current distribution remain largely suitable, while many areas in Central and Northern Europe become less suitable than under the RCP4.5 scenario (Fig. 5c, d). Notably, differences between GCMs are much more pronounced and variation among algorithms is high, particularly under the MIROC GCM (Fig. S3). As under the MIROC RCP8.5 scenario, the GLM model even shows a nearly complete loss of suitable area, we largely discard these results from the discussion.

Discussion
The recent range expansion of the European Bee-eater in Central Europe is often associated with changing climate. In this study, we performed SDMs to test this relationship. Our results confirm a close relationship between climate and realized distribution of the species particularly for the German population. Future predictions point to range extensions at higher latitudes while results also depend on selected model conditions.

Comparison of model approaches
We used eight modelling algorithms to predict the potential current and future distribution of European Bee-eaters.
While we do not recommend to infer generalizations from the results of our single species study in a limited geographical area, our comparison of different model algorithms seems appropriate. Differences in statistical parameters, such as the often criticized AUC values (Lobo et al. 2008;Jiménez-Valverde 2012), were mostly small and output maps had many similarities. Hence, at first glance, an ensemble approach might have been unnecessary. However, we underline the robustness of our results by testing various approaches as has been shown in previous studies (Qiao et al. 2015;Scales et al. 2016). Moreover, a few output maps stand clearly out, in particular, GBM and Maxent, which hardly predict higher latitude populations under current climate. One possible explanation for their deviation in our study might be their different assessment of variable contributions. While mean temperature was the highest contributing thermal variable in all other algorithms, GBM and MAXENT assigned high importance exclusively to minimum temperature (Table 1).
Favouring extreme values over averages may thus lead to more restrictive model outputs: in our case, the predicted unsuitability within the realized distribution of the European Bee-eater in higher latitudes such as Central Europe or the Baltic states. Model outputs vary more distinctly for future predictions, as has been shown elsewhere (Pearson et al. 2006;Ashraf et al. 2017). For instance, the MIROC-GLM-RCP8.5 model predicts a nearly complete loss of suitable areas over the whole breeding range. Again this calls for caution if relying on single model approaches, but offers the opportunity to identify algorithms, GCMs and RCPs that perform less well Thuiller et al. 2019). Each model depends on many factors and our results do not allow to single out a best model approach (Marmion et al. 2009). By performing an ensemble approach, we account for the variation in model outputs which is the general idea behind this methodology (Araújo and New 2007;Elith and Graham 2009;Jones-Farrand et al. 2011;Quillfeldt et al. 2017). A stack of different model outputs has been used for modelling the distribution of a plethora of bird species (e.g. Coetzee et al. 2009;Barbet-Massin et al. 2010;Strubbe et al. 2015;Scales et al. 2016), and we confirm the usefulness of this best practice approach, but also point to model-specific features.

Current potential distribution
The modelled and the realized distribution of the European Bee-eater in the Western Palaearctic are highly similar. Parts of eastern and northern Germany, which were not predicted by all models, have only been colonized by larger numbers of European Bee-eaters within the time frame considered in this study. Hence, the climate data used here could have been inappropriate to describe these fast range dynamics under Fig. 4 Relationship between the area predicted suitable and the population size for a total of 20 European breeding countries and eight algorithms. Country names are given as ISO3166 alpha-2 codes with Germany being highlighted in red. Shaded areas around the regression line highlight 95% confidence limits. Both axes are log (x + 1) scaled rapid global warming. In Central Europe, breeding performance of European Bee-eaters strongly depends on weather conditions (Arbeiter et al. 2016) offering a role for shortterm effects that could conceal effects of climate (Reside et al. 2010). Alternatively, some model approaches, notably Maxent and GBM, underestimated the potential climate niche of the study species. These models partially result in patterns similar to the early assumptions by Voous (1960) who set the northern range edge at the 21 °C July isothermal line with some breeding records below these values. Yet, most algorithms revealed some areas of niche unfilling, i.e. geographic areas that are climatically suitable but uninhabited by the species. These areas are often concentrated along the coasts of Western Europe and the Baltic Sea. Interestingly, there are breeding records from some of these regions (e.g. Birdlife International 2015), but large and persistent colonies are still to be established by the species. While dispersal limitation is often regarded as one range delimiting factor in organisms (Soberón and Peterson 2005), the high mobility of a long-distance migratory bird species known for being a vagrant across many non-breeding areas, (e.g. www. tarsi ger. com, last access 20.01.2020) renders explanation based on dispersal unlikely. However, poleward range shifts in such long-distance migrants are associated with greater migration distances that could also limit further northward range extensions (Doswald et al. 2009;Fiedler et al. 2016;Zurell et al. 2018;Curley et al. 2019). While we acknowledge the great importance of considering the full annual cycle of a migratory species (e.g. Barbet-Massin et al. 2009;Williams et al. 2017;Zurell et al. 2018), data on the wintering behaviour of our study species are far less available than on the breeding distribution, making their consideration a task for future studies. Alternatively, habitat characteristics and reduced food availability, e.g. due to a severe decline of flying insects (Hallmann et al. 2017), could explain some areas of niche unfilling. However, range expansions are a dynamic process that typically occurs uneven over large spatial scales (Veech et al. 2011). Moreover, density-dependent reproductive success known as Allee effects might influence the expansion of a colonial breeding bird species (Serrano et al 2005) but are beyond the scope of our study. Finally, as for all range shifting species, the range expansion of the European Bee-eater conflicts the equilibrium assumptions of SDMs (see Elith et al. 2010 for an overview). SDMs assume that the distribution of a species is in balance with the climatic conditions shaping it. When a species' range shifts, Fig. 5 Ensemble forecasts for the year 2050 using two global circulation models (CCSM and MIROC) and two relative concentration pathways (RCP 4.5 and 8.5). Shading represents the number of algorithms predicting potential occurrence within the background area either favourable climatic conditions stretch far beyond what a species was able to colonize yet (i.e. creating a time lag in species response to climate change) or a species expands into novel climatic combinations not reflected in the original range but which are still favourable (e.g. as for many invasive species). This disequilibrium in the species-environment relationship in the European Bee-eater may explain underestimated regions of the realized distribution but also the strong variation in future predictions (see below, Elith et al. 2010). In sum, no matter whether the source of uncertainty in the predictions arises for methodological or biological reasons, overall model outputs still underline climate as main driver of the European and German current distribution of the European Bee-eater.

Relationship between climate suitability and population size
The positive relationship between climate suitability and population size among European countries is in concordance with our hypotheses that the breeding distribution of the European Bee-eater is mainly defined by climatic factors. This includes the German population which fits within or is close to the 95% confidence interval. Besides this general relationship, however, we found that confidence limits stretch widely-often by two orders of magnitude in most model approaches. This is most likely due to a combination of uncertainties in the predicted suitable areas from the respective SDM but also in the country-level population estimates. It is noteworthy that some model approaches predicted population size appropriately, but were not able to reveal the underlying spatial pattern. For instance, FDA, GBM, GLM, and Maxent failed to predict the large populations in eastern Germany and potentially overestimated population sizes in southwestern Germany. ANN and RF show patterns of overprediction as many areas defined as suitable contrast to the realized distribution of European Bee-eaters in Germany. From this, while the amount of suitable area in Germany is in line with the population estimate, the location of suitable areas often diverges from the realized distribution of European bee-eaters in Germany calling for a closer inspection of the colonization history in the country and using climate information of finer temporal resolution. On a European scale, the positive correlation between population size and climate suitability for Belarus and the Baltic states of Latvia and Lithuania was absent for some SDMs. This is possibly due to a lower data coverage in Eastern Europe and Asia compared to Western and Southern Europe (Fourcade et al. 2014b). Belarus was in almost any algorithm far outside the confidence limits (Fig. 4) and excluding this single country alone from the country-level analysis of the amount of suitable area from population size improved that relationship tremendously (data not shown).
Besides methodological effects of data coverage and quality, spatial variation of predictor variable collinearity could be a further reason that would call for further investigation. However, we confirm the general positive relationship of environmental factors and abundance (VanDerWal et al. 2009;Thuiller 2014) and add to the studies that infer data on avian population size and density from SDMs (Oliver et al. 2012;Barker et al. 2014;Carrascal et al. 2015).
Recent declines in southern parts of the range as observed e.g. in Turkey (Birdlife International 2015) occur in areas that show suitable conditions under current and most future projections. Hence, we hesitate to link these patterns to direct effects of global warming. Instead, other factors like agricultural intensification and bird hunting should not be neglected as possible explanations (Donald et al. 2001;Brochet et al. 2016). Moreover, European Bee-eaters need large flying insects as food which have shown marked declines in recent years and are particularly affected by climate change (Hallmann et al. 2017;Soroye et al. 2020) offering the potential for indirect negative effects on populations of European Bee-eaters.

Future potential distribution
Future projections under the RCP4.5 emission scenario predict a shift of the northern range edge as has already been shown for many species (e.g. Thomas and Lennon 1999). Areas in Sweden, where the species already bred in the past (Cramp 1985), and southern Finland will have suitable climate, while predictions for large parts of the United Kingdom, where the species already bred (e.g. Holling and Rare Breeding Birds Panel 2017), vary depending on selected algorithms. Future projections under the RCP8.5 emission scenario showed a large variation in suitable climate according to the two GCMs considered. However, while the CCSM RCP8.5 showed a level of consistency among SDM algorithms compared to the RCP4.5 GCMs, the MIROC RCP8.5 predictions deviated strongly from each other-up to a point where the entire current distribution was predicted unsuitable (in the GLM, Fig. S3). One possible explanation is that climatic changes diverge strongest among GCMs in RCP8.5 as this represents the worst-case scenario in future global climate developments (Hausfather and Peters 2020). Under these extreme conditions, some GCMs could predict locally much stronger changes than others-with profound effect of SDM projections. Differences between circulation models have also been shown for other projections, stressing the need for considering more than one GCM (Schidelko et al. 2011(Schidelko et al. , 2013Thuiller et al. 2019). However, further research is needed to compare how strong and in which geographic areas and under which RCP GCMs are diverging most. In contrast to other bird species, losses of suitable breeding range of the European Bee-Eater were not concentrated on the southern range edge (Thomas and Lennon 1999;Virkkala and Lehikoinen 2014), and these areas remain climatically suitable under most predictions. Based on climate projections, our study species could well survive in the Maghreb area, which has been previously shown to be important for future avian diversity under climate change (Barbet-Massin et al. 2010). Hence, our hypothesis that global warming will heavily affect the species in the southernmost parts of its range is not supported by our results. European Bee-eaters might be tolerant to predicted rising temperatures, perhaps due to plasticity in thermoregulatory behaviour (Glutz von Blotzheim 1980;Yosef 2010). While the species even breeds in desert areas, it is typically confined to river valleys and other specific habitats, which were not considered within our approach.

Conclusions
Taken together, the European Bee-eater might indeed be an iconic messenger of climate change. At least in the short-term and under medium emission scenarios, European Bee-eaters are predicted to shift northwards, while currently populated areas remain suitable, especially in the south of the distribution. More regionally, the current population size in Germany can be predicted by climatic conditions during the breeding season by most algorithms and is in line with conditions in other European countries. Therefore, the European Bee-eater can be regarded as an appropriate 'indicator species' for global warming in Central Europe (see Gregory et al. 2009), and we emphasize the need for species-specific analyses of European species' distributions under climate change. Future studies on regional short-term population dynamics of the European Bee-eater should consider weather (Arbeiter et al. 2016) and the availability of suitable breeding habitats.
Author contributions All authors contributed to the ideas of the project. Analysis was led by JOE with critical input from DS. DS led the writing with support from JOE and KS. All authors contributed to the drafts and gave final approval for publications. The authors declare no conflict of interest.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.