Global warming promotes biological invasion of a honey bee pest.

Climate change and biological invasions are two major global environmental challenges. Both may interact, e.g. via altered impact and distribution of invasive alien species. Even though invasive species play a key role for compromising the health of honey bees, the impact of climate change on the severity of such species is still unknown. The small hive beetle (SHB, Aethina tumida, Murray) is a parasite of honey bee colonies. It is endemic to sub-Saharan Africa and has established populations on all continents except Antarctica. Since SHBs pupate in soil, pupation performance is governed foremost by two abiotic factors, soil temperature and moisture, which will be affected by climate change. Here, we investigated SHB invasion risk globally under current and future climate scenarios. We modelled survival and development time during pupation (=pupal performance) in response to soil temperature and soil moisture using published and novel experimental data. Presence data on SHB distribution were used for model validation. We then linked the model with global soil data in order to classify areas (resolution: 10 arcmin; i.e. 18.6 km at the equator) as unsuitable, marginal and suitable for SHB pupation performance. Under the current climate, the results show that many areas globally yet uninvaded are actually suitable, suggesting considerable SHB invasion risk. Future scenarios of global warming project a vehement increase in climatic suitability for SHB and corresponding potential for invasion, especially in the temperate regions of the Northern hemisphere, thereby creating demand for enhanced and adapted mitigation and management. Our analysis shows, for the first time, effects of global warming on a honey bee pest and will help areas at risk to prepare adequately. In conclusion, this is a clear case for global warming promoting biological invasion of a pest species with severe potential to harm important pollinator species globally.

INTRODUCTION

Human‐mediated biological invasion is considered to be one of the most serious threats for biodiversity (Dyer et al., 2017; Early et al., 2016; McGeoch et al., 2010), which may even cause the breakdown of classical biogeographic regions (Capinha, Essl, Seebens, Moser, & Pereira, 2015). Moreover, invasive species can cause considerable social, economic and ecological damage (Blackburn et al., 2011; Pimentel, Lach, Zuniga, & Morrison, 2000), altering ecosystems and endangering food security (Schweiger et al., 2010; Veldtman et al., 2011; Ziska, Blumenthal, Runion, Hunt, & Diaz‐Soltero, 2011). The invasiveness and impact of invasive species is a complex interplay between biotic and abiotic factors (D'Antonio, 1993; Thuiller, Richardson, Rouget, Procheş, & Wilson, 2006; Tobin, 2015) with varying consequences (Cuthbert, Dickey, McMorrow, Laverty, & Dick, 2018; Rejmánek & Richardson, 1996; Ricciardi & Cohen, 2007). Changing climates can impact the current status of alien species, often resulting in an increased probability to become established or to spread to areas currently deemed environmentally unsuitable (Dukes & Mooney, 1999; Early et al., 2016; Sutherst, Floyd, & Maywald, 1996). This is particularly true for ectotherms, which depend on climatic conditions to permit survival and development within the thermal limitations a habitat poses (Barbet‐Massin et al., 2013; McCann, Greenlees, & Shine, 2017; Roura‐Pascual et al., 2004).

The impact of invasive species is of major concern for society when the provision of ecosystem services is affected (Pejchar & Mooney, 2009). The Western honey bee, Apis mellifera, Linnaeus, is a particularly important species for providing pollination services globally (Calderone, 2012; Hung, Kingston, Albrecht, Holway, & Kohn, 2018). However, managed honey bees have been facing severe colony losses in recent decades (Brodschneider et al., 2018; Jacques et al., 2017; Neumann & Carreck, 2010; Van Engelsdorp & meixner, 2010). Even though the number of managed honey bee colonies is increasing globally, the demand for pollination is growing at a much higher rate (Aizen & Harder, 2009; Gallai, Salles, Settele, & Vaissière, 2009). Amongst the many factors potentially impacting honey bee health and thus pollination services, invasive parasitic species, e.g. introduced by global trade in honey bees and related products such as wax and honey (Chanpanitkitchote et al., 2018; Krongdang, Evans, Chen, Mookhploy, & Chantawannakul, 2018; Neumann, Pettis, & Schäfer, 2016; Ouessou Idrissou, Huang, Yañez, & Neumann, 2019; Schäfer et al., 2019), can play a key role (Neumann et al., 2016; Potts, Biemeijer, et al., 2010; Rosenkranz, Aumeier, & Ziegelmann, 2010; van Dooremalen, Cornelissen, Poleij‐Hok‐Ahin, & Blacquière, 2018). However, knowledge of the potential effects of climate change on such species is currently lacking (Le Conte & Navajas, 2008). The small hive beetle (SHB, Aethina tumida, Murray) is a long known parasite of social bee colonies (honey bees: A. mellifera (Lundie, 1940); Apis cerana (Cervancia, Guzman, Polintan, Dupo, & Locsin, 2016), bumblebees: Bombus impatiens (Spiewok & Neumann, 2006), stingless bees: (Greco et al., 2010) native to sub‐Saharan Africa, which can also infest nests of solitary bees (Megachile rotundata, Gonthier et al., 2019). Since 1996, SHB has become an invasive species and has established local populations on every continent except Antarctica (Neumann et al., 2016; Schäfer et al., 2019). Despite comprehensive elimination and contingency efforts, it is likely to continue spreading (Schäfer et al., 2019). The impact of SHBs on honey bee colonies in the invasive ranges is well documented (Neumann & Elzen, 2004) and depends on infestation levels, with higher infestation levels more likely leading to host colony collapse (Spiewok et al., 2007). When SHBs mass reproduce, with often thousands of larvae (Neumann & Elzen, 2004), they can kill even strong colonies of European honey bee subspecies within 10 days (Neumann, Hoffmann, Duncan, & Spooner‐Hart, 2010), often resulting in the full structural collapse of the entire nest (Hepburn & Radloff, 1998). This is very rare in the native range of SHB in Africa in colonies of the respective local honey bee subspecies (Lundie, 1940; Neumann, 2017; Schmolke, 1974), where SHBs probably mostly rely on non‐destructive low‐level reproduction (Ouessou Idrissou, Straub, & Neumann, 2018). The higher susceptibility of European honey bee subspecies is probably due to quantitative differences in a range of social immunity traits compared to the African ones (e.g. aggression; Elzen et al., 2001, absconding [non‐reproductive swarming; Neumann et al., 2018] and social encapsulation of SHBs; Neumann et al., 2001).

Besides biotic factors, abiotic factors may also contribute to the invasion success of SHBs. In contrast to other beetles, which can complete an entire life cycle within host colonies (Krishnan, Neumann, Ahmad, & Pimid, 2015), SHBs have to pupate in the soil to complete their life cycle (Ellis, Hepburn, Luckman, & Elzen, 2004; Lundie, 1940). SHB pupation success (survival rate) and the duration of pupation are governed by soil humidity and temperature (Akinwande & Neumann, 2018; Bernier, Fournier, & Giovenazzo, 2014; Ellis et al., 2004; Meikle & Diaz, 2012; Meikle & Patt, 2011). It is therefore apparent that abiotic factors can play a key role in explaining the performance and thus invasion success of this species. Indeed, under favourable environmental conditions, i.e. high humidity and temperature, SHBs can cause significant damage to apiculture outside its endemic range. For instance, in 1998 SHB caused damage of more than 3 million USD in Florida (Neumann & Elzen, 2004). There may be up to six SHB generations per year under US and South African climatic conditions (Neumann & Elzen, 2004), and De Guzman and Frake (2007) showed that almost 16 complete life cycles can be achieved within a year under a constant soil temperature of 34°C.

Given the importance of social bees for pollination services and their economic value (Gallai et al., 2009; Klein et al., 2006; Velthuis & Van Doorn, 2006), assessing the risks of SHB invading currently uninvaded areas and potential changes in the severity of SHB impacts is of utmost importance, especially under changing climatic conditions. This information is urgently needed (EFSA, 2015), as it will define management strategies during different stages of invasion (Cook, Thomas, Cunningham, Anderson, & Barro, 2007; Schäfer et al., 2019). Therefore, identifying environmental limitations and their changes are key to assess the invasiveness of alien species and their biotic interactions (Schweiger et al., 2010).

Applying a common correlational approach of assessing the climatic niche of a species based on distributional data (Araújo & New, 2007; Thuiller, Lafourcade, Engler, & Araújo, 2009) is particularly difficult for alien species that are still spreading and thus not in equilibrium with the environment and for which data from the native range are scarce (Václavík & Meentemeyer, 2012). Since this all applies to SHB, we developed a mechanistic niche model relying on physiological tolerances to environmental conditions and the corresponding effects on performance. However, the impact of environmental factors can vary among life stages, and thus it is important to focus on the most sensitive ones (Bowler & Terblanche, 2008). The part most sensitive to environmental conditions during the life cycle of SHB is the pupal stage outside the host colony and therefore we focused on measures of pupal performance. We used empirical data on the response of survival rate and development time to soil moisture and temperature conditions and assessed the global invasiveness and severity of SHB under current and projected future climatic conditions. We assume that pupal performance is one key aspect related to the invasion risk of SHB at a global scale and we predict that SHB invasion risk will increase as climate change and global warming in particular, promotes the chances of SHB to survive and thrive in many areas of the world.

MATERIAL AND METHODS

SHB pupal performance data

Small hive beetle pupal performance data are here defined as survival rate and developmental time and were collected from peer‐reviewed literature with focus on the impact of soil temperature and moisture (Table S1). Since soil moisture was provided either as weight (moisture) or volume ratios (soil water content), the gravimetric measures were converted into volumetric measures according to the bulk density (kg/m3) of the used soil types (Table S1). To fill identified data gaps in the published studies, additional laboratory experiments were performed (Table S1 and Method S1).

SHB pupal performance curves

To assess the potential global distribution and invasiveness under given and projected future climatic conditions, we quantified the responses of pupal survival rate and developmental time to varying soil temperatures and moisture conditions and combined them into a composite measure of pupal performance. Performance curves for many physiological processes are well described and they typically rise to an optimum and then decline more or less steeply to zero performance (e.g. Huey & Kingsolver, 1993). For survival rate and soil temperature, we applied a performance function (S(T)), used for ectothermic invertebrates, e.g. in Deutsch et al. (2008) or Vasseur et al. (2014), where the rise is described by a Gaussian function and the decline by a parabolic function:where S is the performance metric of pupal survival, T is the soil temperature, T opt is the temperature with maximum performance, T max is the upper critical temperature at which performance is zero and σ p is a shape parameter determining the steepness of the Gaussian function.

For survival rate and soil moisture (S(M)), we used the same function, but due to a larger plateau in the response curve, we had to re‐parameterise the exponent of the parabolic function. Therefore, we let the exponent increase from 2 onwards and selected the best fitting model according to the lowest value of Akaike's information criterion (AIC; see Figure S1 for visual assessment and Table S2 for AIC values) leading to the following model:where S is the performance metric of pupal survival, M is the soil moisture, M opt is the moisture with maximum performance, M max is the upper critical moisture at which performance is zero and σ p is a shape parameter determining the steepness of the Gaussian function.

To allow for varying shapes of the thermal performance curve under different soil moisture conditions and vice versa, we combined both performance curves (S(TM)) via an interaction term:where z is a scaling factor.

We used a nonlinear regression approach to derive performance curves by fitting observed pupal survival rates (S(T) and S(M)) at the respective temperature (T) and moisture conditions (M) to the Equations (1.1) and (1.2) and estimated the parameters T opt, T max, M opt, M max and σ p. Starting values for the iterative estimation approach, by minimizing sum of squares, were obtained by visual inspections of plotting survival rate against temperature and soil moisture. We compared models including only one function (either S(T) or S(M)), models including their additive and their interactive effects. The lowest AIC values indicated that the interactive effect performed best (Table S2).

Development times often follow a u‐shaped relationship with temperature, but if there is no indication of an increase in development times at very high temperatures, an asymptotic exponential function can be used (see e.g. Kingsolver, Diamond, & Buckley, 2013). We used a three‐parameter asymptotic exponential function (D(T)) for both soil temperature and soil moisture:where a is the horizontal asymptote on the right hand side (i.e. at high temperatures) and defines the minimum number of days for pupal development at high temperatures, b is given by a – R0, where R0 is the intercept, and c is the rate constant defining the shape of the curve. Parameters were also estimated with nonlinear regressions. Model comparison based on AIC indicated that only soil temperature is relevant for development time (Table S2).

Finally, we combined performance measures of pupal survival rate and development time into a composite measure of pupal performance (P) by:where dividing by the maximum we let the measure vary between 0 and 1.

Global soil temperature and moisture data

We predicted global distribution and invasiveness of SHB under current climatic conditions using the composite measure of pupal performance (P) based on environmental information on soil temperature and soil moisture. For soil temperature, we used data provided by microclim (Kearney, Isaac, & Porter, 2014). Microclim provides hourly estimates from the surface to 1 m depth for the middle day of each month at a resolution of 10 arcmin (i.e. 18.6 km at the equator) including six shade levels and three substrate types (soil, rock and sand). According to the biology of SHB (De Guzman, Frake, & Rinderer, 2010; Pettis & Shimanuki, 2000), we extracted data for soil at a depth of 10 cm and for a subset of every third hour. To assess the level of shading by vegetation cover, we used data on the Normalized Difference Vegetation Index (NDVI) as a proxy. NDVI data were obtained from the Global Inventory Modeling and Mapping Studies (GIMMS; Pizon, 2005; Tucker et al., 2005) of the years 1981–2010 (Tucker, Pizon, & Brown, 2016) at a biweekly interval and at a resolution of 5 arcmin (9.3 km at the equator). These data were aggregated to mean monthly values at the 10 arcmin resolution of the soil temperature data and averaged across the 30‐year period. NDVI values were equally binned into five classes and assigned to the respective levels of shading for the microclim data (0%, 25%, 50%, 75%, 100%). Based on the NDVI values, we extracted the respective soil temperature data for each grid cell leading to global soil temperature estimates for the middle day of each month for every third hour.

Soil moisture data were obtained from the ESA CCI Surface Soil Moisture (ESA CCI SM) project v2.2 (Liu et al., 2011, 2012; Wagner et al., 2012) on a daily basis for the period from 1985–2014 at a resolution of 15 arcmin (27.9 km at the equator). We calculated monthly means across the 30 year period and disaggregated the data to the 10 arcmin resolution. Since ESA CCI SM data did not cover the tropics, we filled these gaps with data from NASA SMAP L4_SM data product (Reichle, Lannoy, Koster, Crow, & Kimball, 2017) on a 3‐hourly basis for the period from 2015–2017 at a 9 km resolution (provided on a global cylindrical equal area grid). We calculated mean monthly values for the 3‐year period and after re‐projection we aggregated the data to the 10 arcmin resolution.

Since SHB has only been observed in vegetated areas we masked non‐vegetated areas using the NASA Land Cover Type Climate Modelling Grid product (MCD12C1; Nasa Lp Daac, 2013). These data on dominant land cover types, originally provided at a 3 arcmin (5.6 km at the equator) resolution, were again aggregated to the 10 arcmin grid.

To assess potential consequences of global warming on the future invasion of SHB, we used current mean monthly surface temperature averaged across the period from 1960–1990 from WorldClim (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005) provided at the 10 arcmin resolution and future projections obtained from the HadGEM2‐ES general circulation model of the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (Stocker, 2014). We used two scenarios of representative concentration pathways (RCPs) for 2060 (averaged over 2041–2060) and 2080 (averaged over 2061–2080) resulting in an average global increase of 2.88°C and 2.98°C (RCP2.6) or 4.24°C and 6.09°C (RCP8.5). To estimate future soil temperatures, we first performed linear regression models for each of the 3‐hourly soil temperature data as a function of mean surface temperature in a particular month. In this way we captured diurnal and seasonal variation in these relationships. We found strong relationships (mean R 2 = 0.97, range = 0.90–0.99) but with varying slopes (mean = 1.00, range = 0.60–1.13) and intercepts (mean = 0.44, range = −2.22–12.21). Subsequently, we used these models to predict future soil temperatures per grid cell for each 3‐hr interval for the middle day of each month. To preserve grid‐cell specific deviations from these overall relationships, e.g. caused by slope, aspect or precipitation, we added the respective residuals from the regression analyses under current conditions to the projections under future conditions, assuming that these deviations are constant over time. We also tested this approach for soil moisture and its relationship to monthly precipitation, but very low R 2 values (mean = 0.13, range = 0.04–0.24) indicated low reliability. We therefore relied on scenarios of temperature change only.

Global predictions and future projections of SHB pupal performance

We used performance curves of pupal survival rates (S(TM)) and development time (D(T)) to predict both processes separately based on current soil temperature and moisture conditions. To account for the strong geographic differences in diurnal variation of soil temperature and the corresponding consequences for pupal performance we predicted performance for each of the 3‐hr intervals per month and integrated them in a second step by averaging. We then used Equation (3.1) to calculate the composite index of pupal performance per month, whereas we used the maximum of S(TM)/D(TM) across all grid cells and months (global maps for monthly survival rate, development time and pupal performance are provided in Figures [Link], [Link], [Link]). Monthly pupal performance was further condensed in two ways: (a) it was averaged across the months per grid cell assuming that pupal performance accumulates across the varying conditions within a year; and (b) by extracting the highest level of performance across the months per grid cell assuming that invasiveness depends on maximum performance during shorter periods.

To assess the predictive ability of the model in general and to discriminate the relevance of mean annual climatic conditions from short‐term optimal conditions for pupal performance, we used actual reported georeferenced occurrences in the native and invaded range (Table S3). We included only established populations by focusing on observations of 3 years or longer. Since these data represent presence‐only data, we used the continuous Boyce index (Boyce, Vernier, Nielsen, & Schmiegelow, 2002; Hirzel, Lay, Helfer, Randin, & Guisan, 2006) to assess the quality of our predictions. This index varies from −1 (worse than expected by chance) to 0 (not better than expected by chance) to 1 (perfect predictions). The Boyce index compares the predicted frequency distribution of evaluation points with their expected frequency based on the distribution within a selected area and is thus sensitive to the spatial extent of the selected area. To overcome a potential bias by selecting a too large area, we analysed occurrences in different regions separately by calculating convex hulls for three areas in Africa, three in North America and one in Australia with sufficient data points. Further, we considered some uncertainty in the georeferences of the observations and used the maximum of pupal performance across all grid cells within a 25 km buffer. For evaluation of the Boyce index, we used an analogy to the categorisation recommended by Landis and Koch (1977) for the true skill statistic, which also ranges between −1 and 1, the following: excellent, Boyce index > 0.75; good, 0.40 < Boyce index < 0.75; and poor, Boyce index < 0.40.

Future projections of pupal survival rates, development time and pupal performance were calculated analogous to current conditions but for means of comparison, we used the same value of max[S(TM)/D(TM)] for Equation (3.1) as for current conditions.

All analyses were performed in the statistical environment R (R Core Team, 2016) using the packages colorRamps (Keitt, 2012), ecospat (Broennimann, Di Cola, & Guisan, 2016), gdalUtils (Greenberg & Mattiuzzi, 2015), gtools (Warnes, Bolker, & Lumley, 2015), maptools (Bivand & Lewin‐Koh, 2016), minpack.lm (Elzhov, Mullen, Spiess, & Bolker, 2016), ncdf4 (Pierce, 2015), raster (Hijmans, 2017), RColorBrewer (Neuwirth, 2014), rgdal (Bivand, Keitt, & Rowlingson, 2016) and sp (Pebesma & Bivand, 2005).

RESULTS

Climatic factors defining pupal performance of SHB

Models for both SHB pupal survival rate and developmental time fitted the experimental data very well (Figure 1). Interestingly, we found an unusual positive skew of the thermal performance curve for survival rate with a very steep slope for the rising part of the curve and a less steep decreasing part (Figure 1a). Estimated temperature optimum for survival rate was at 18.9°C and the upper critical temperature was 40.9°C (all parameter estimates are provided in Table S4). Pupal survival rate of SHB also showed a broad range of suitable soil moisture conditions reaching optimal conditions already at 0.04 m3/m3, followed by a larger plateau and reaching critical moisture conditions at 0.77 m3/m3 (Figure 1b). Pupal development time was well described by the asymptotic exponential relationship with temperature, with the asymptote approaching a minimum time for development of 15 days (parameter a in Table S4). Survival rate and development time indicated a temperature‐related trade‐off with longer development times at optimal temperatures for survival. The composite index of pupal performance consequently led to a shift of optimal temperature from 18.9°C for survival rate to 27.5°C when development time was additionally considered (Figure 1d). At this temperature, development time decreased from 43 to 18 days (Figure 1c), while survival rate was still 87% (Figure 1a).

Figure 1

Performance curves of small hive beetle pupal survival rate in response to temperature (a) and soil moisture (b), pupal development time (c) and composite pupal performance (d) as a function of temperature. Since pupal survival was calculated based on interactive effects between thermal and moisture performance (Equation 1.3), visualised curves are based on constant average conditions of the respective environmental measure, i.e. soil moisture of 0.21 mm3/mm3 for thermal performance (a) and soil temperature of 25.8°C for moisture performance (b) and respective upper and lower 10% quantiles of the data points are not displayed, i.e. data points with values lower than 0.02 and values higher than 0.56 mm3/mm3 for thermal performance (a) and data points with values lower than 16.0°C and values higher than 37.6°C for moisture performance (b). Black dashed lines in (a) and (b) indicate optimal conditions (T opt and M opt in Equations 1.1. and 1.2), red dotted lines indicate the upper critical temperature T max (a) and the upper critical moisture M max (b). Green dashes lines indicate the optimal temperature for composite pupal performance (d) and visualise a decrease in development time from 43 days according to T opt to 18 days (c) while survival rate decreased only to 87% (a). Point colours indicate data sources: black, Bernier et al. (2014); red, Ellis et al. (2004); green, Meikle and Diaz (2012); dark blue, Meikle and Patt (2011); light blue, unpublished data from new experiment 1; magenta, unpublished data from new experiment 2 (see Table S1 and Method S1)

Boyce indices for the selected regions and all regions combined show that highly beneficial, but short‐term climatic conditions, i.e. conditions in the ‘best’ month, explain the distribution and invasiveness of SHB much better than long‐term conditions averaged across the year (Table 1). Predictive ability was generally good with one exception for the region around Lake Victoria in Africa.

Table 1

Boyce index calculated for three regions in Africa (af), three regions in North America (na), one region in Australia (aust) and all regions combined (all; see Figure 3)

Region	Average	Best
af1	−0.44	0.67
af2	0.58	0.66
af3	−0.55	−0.72
na1	0.74	0.83
na2	0.67	0.46
na3	1.00	1.00
aust	0.94	0.76
all	−0.02	0.70

Average: Boyce index computed for pupal performance measures based on annual averages, assuming that overall pupal performance accumulates across the varying conditions within a year.

Best: Boyce index computed for pupal performance measures based on the highest value across the months, assuming that invasiveness depends on maximum performance during shorter periods.

Figure 3

Global predicted pupal performance of small hive beetle (SHB) for Europe (a), Asia (b), Africa (c), Australia and the Malay Archipelago (d), North and Central America (e) and South America (f). Pupal performance is based on a composite index combing pupal survival rate and development time (Equation 3.1) and ranges between zero (no performance) and one (maximum performance). According to thresholds obtained from model validation (see Figure 2), continuous pupal performance values were classified into conditions of high climatic suitability (values higher than 0.64; red to orange colours), marginally suitable (values between 0.4 and 0.64; yellow to green) and unsuitable climatic conditions (values below 0.4; blue to grey colours). Non‐vegetated areas are masked in white. Open circles show locations with georeferenced occurrences of SHB. Black polygons depict areas used to determine expected frequency distribution of performance values used to assess predictive ability with the Boyce index (see Table 1)

For a better interpretation of the continuous pupal performance index, we identified two thresholds using the predicted/expected ratio used to calculate the Boyce index (Figure 2). A predicted/expected ratio higher than one indicates a better prediction than expected by chance (Hirzel et al., 2006) and occurred at pupal performance values higher than 0.64. In such areas climatic conditions are considered as highly suitable.

Figure 2

Predictive ability of the small hive beetle pupal performance model. Predictive ability is measured by the ratio of the frequency of actually predicted performance values for the evaluation points over the expected frequency based on the predictions for the entire considered area, in this case for all evaluation areas (see Figure 3) combined (see Hirzel et al., 2006). Horizontal grey dotted line indicates a ratio of one. Above this line, pupal performance is predicted more often than expected by chance. Areas with performance values above the corresponding threshold of 0.64 (right dashed vertical line) are considered as highly suitable for pupation. Below a threshold of 0.4 (left dashed vertical line) no established populations were reported and areas with such predicted performance values are considered as unsuitable. The intermediate range is considered as marginally suitable

The observed drop in the predicted/expected ratio at predicted pupal performance higher than 0.9 might be caused by a systematic bias in the observations (too few observations in the “best” areas) or caused by other environmental or biotic conditions apart from climate. We further observed SHB at predicted pupal performance values between 0.64 and 0.40, but at lower frequencies than expected by chance. Conditions in such areas are considered as marginal. No observations have been made in areas of predicted pupal performance values lower than 0.40 and are thus considered as unsuitable.

Predictions of pupal performance of SHB under current climatic conditions indicated high climatic suitability in its native range in sub‐Saharan Africa and generally in the Southern hemisphere (Figure 3). Here, all (sub‐)alpine areas, almost all of New Zealand and the southern‐most part of South America, below 45° South latitude, are unsuitable for pupation (Figure 3). Climatic constraints currently limit the distribution of SHB for large parts of the Northern hemisphere. Marginal to optimal conditions for pupation occur in vegetated areas in Asia, North Africa, Southern Europe and North and Central America up to 57° North latitude, with some exceptions up to 60° North latitude (Figure 3).

Under a moderate warming scenario (RCP2.6), large areas in the Northern hemisphere are projected to become highly suitable by 2060 (28.5% gain; Supplement Table S5) but remaining rather constant until 2080 (28.8% gain; Figure 4a,b). The northern boundary for marginal pupation performance is projected to shift to 67°N in 2080. In particular, large areas in the Russian Federation, Canada and Europe could become suitable for SHBs to pupate and thus establish populations. Moreover, in areas where SHB pupal performance under current conditions could be considered marginal, conditions are likely to shift to optimal. No major changes in pupation performance are projected for the southern hemisphere (Table S5), with a slight boundary shift southward from 45°S to 46°S.

Figure 4

Pupal performance of small hive beetle projected to the representative concentration pathways (RCPs) 2.6 (a,b) and 8.5 (c,d) for the years 2060 (a,c) and 2080 (b,d). Pupal performance is based on a composite index combing pupal survival rate and development time (Equation 3.1) and ranges between zero (no performance) and one (maximum performance). According to thresholds obtained from model validation (see Figure 2), continuous pupal performance values were classified into conditions of high climatic suitability (values higher than 0.64; red to orange colours), marginally suitable (values between 0.4 and 0.64; yellow to green) and unsuitable climatic conditions (values below 0.4; blue to grey colours). Non‐vegetated areas are masked in white

Under the more severe warming scenario (RCP8.5), projected increases in pupal performance are drastic for the Northern hemisphere (2060: 48.2% gain, 2080: 84% gain; Figure 4c,d). By 2060, the northern boundary for marginal pupal performance shifted to 69°N and to 71°N by 2080, projecting the possibility of marginal SHB pupal performance on the South Island of the Novaya Zemlya archipelago. In the Southern hemisphere, changes are still less pronounced (Table S5). A maximum southward shift of marginal to optimal performance to 49°S is projected for 2080. Furthermore, large areas of New Zealand and Tasmania will become suitable for SHB pupation.

DISCUSSION

Here we present the first study to assess the impact of global warming on an invasive honey bee pest on a global scale. Our results show a high invasion risk across the globe with potential dire consequences for its hosts. Moreover, the risk increases considerably with increasing temperatures in the future. By categorising pupal performance from unsuitable to marginal to optimal, we show that SHB can potentially colonize an area much larger than is currently the case, confirming earlier concerns (Neumann et al., 2016).

With our mechanistic approach, combining impacts of soil temperature and moisture on pupal survival and development time into a composite thermal performance curve for pupation, we can go beyond a mere expression of survival and define its importance in relation to habitat suitability (Kearney & Porter, 2009) and the impact on bees. For instance, development time is a limiting factor for SHB performance under temperate climatic conditions, since short growing seasons can prevent the completion of metamorphosis (Bernier et al., 2014). Under warmer climatic conditions, developmental time is not a limiting factor for pupal survival, but it can foremost be a predictor of the number of life cycles that can be completed during a growing season, which is indicative of the population build‐up and moreover of the consequences for bees (Ellis et al., 2004; Neumann et al., 2016; Spiewok et al., 2007).

Performance curves of ectotherms, most often measured as thermal response curves, generally show a slowly increasing slope to optimal performance and a sharp decrease thereafter (Dowd, King, & Denny, 2015; Kingsolver et al., 2011), partially due to structural constraints (Huey & Kingsolver, 1989). The performance curve of SHB pupal survival could be considered an unusual shark fin shape, with a sharp increase in survival from the lower thermal limit over a short temperature interval and a gradual decrease in survival towards the upper thermal limit. A similar survival curve has been shown for other beetle species pupating in the soil (Entomoscelis americana, Lamb & Gerber, 1985) and other ectotherms like amphibians (Bachmann, 1969). Our results indicate that SHB is highly sensitive to small changes at the lower end of its temperature niche, which is similar to other beetle species pupating in soil such as E. americana pupae exposed to constant low temperatures in the laboratory (Lamb & Gerber, 1985). However, under field conditions temperatures are fluctuating and may lead to a different result. For instance, E. americana larvae and pupae were able to survive lower temperatures and also used thermo‐regulation to partially overcome fluctuations and constraining temperatures (Lamb & Gerber, 1985). As a consequence, insect species may survive in cooler temperatures under natural conditions than explained by their thermal niche modelled on the basis of constant laboratory conditions. Further, temperatures can also spatially vary within the grid cells we used for modelling. While our models predict general patterns of pupal performance at a global scale, specific local microclimatic conditions, e.g. heat island effects of cities in cooler areas, might, in some cases, allow for sustainable local populations even in areas currently marked as unsuitable. However, the reproductive potential of such populations might be significantly reduced because of the limited time available for completing one or more generation cycles.

Overall, predicted SHB pupal performance fitted very well with observed presence data of SHB. Occurrences in New South Wales, Australia and Florida, USA and even latitudinal outliers like Ontario, Canada are well explained (see Figure 3d,e; Table 1). Nevertheless, invasions beyond the currently predicted limits of SHBs distribution might also be possible if the species manages to adapt to novel conditions (Atwater et al., 2016; Chapman, Scalone, Štefanić, & Bullock, 2017; Krehenwinkel, Rödder, & Tautz, 2015). Since SHBs naturally occur from the Kalahari to equatorial rainforests of sub‐Saharan Africa (Neumann & Elzen, 2004), it is inevitable that the native range holds different ecotypes (Neumann et al., 2016) likely leading to a high adaptive potential of this species after invasions.

Combining models of pupal survival and development time led to a shift of optimal temperature from 19°C for survival to 28°C for overall pupal performance. This shift is predictable given that only a marginal difference in the survival rate between the thermal optimum (19°C) and high temperatures up to 30°C can be observed while developmental time exponentially decreases with increasing temperatures (Figure 1). The strong increase of SHB pupal survival rate in response to small changes at lower temperatures (Figure 1a) is also reflected by a high sensitivity of overall pupal performance (Figure 1d). This high sensitivity is the likely cause of the considerable increases of climatically suitable areas even under a moderate warming scenario (RCP2.6). The rather tropical nature of the thermal niche profile, compared to a temperate one (Deutsch et al., 2008; Huey et al., 2012), on the other hand, can explain why no areas were projected to become unsuitable towards the upper limit of the thermal tolerance of SHB even under the most severe scenario. Furthermore, areas where currently only marginal pupal performance is predicted are likely to facilitate optimal performance with increasing temperatures. For future projections of SHB pupal performance, we could not include scenarios of changing soil moisture (see Section 2), but since SHB has a very broad tolerance level for soil moisture conditions and is extremely insensitive to dry soils (Table S1), we believe that projected soil moisture anomalies would not have impacted our general results appreciably.

Small hive beetle can occur in particularly dry, sparsely vegetated semi‐desert areas such as the Kalahari in Namibia and Botswana when hosts are present (Ellis & Munn, 2005; Phokedi, 1985). Nonetheless, such extreme environments pose a challenge for beetles to survive (El‐Niweiri, El‐Sarrag, & Neumann, 2008). In particular soil physics might prevent SHB to successfully complete the life cycle. Under laboratory conditions, SHB has been shown to successfully pupate below 3 cm (Meikle & Diaz, 2012) and it has been found at depths up to 20 cm in the field (Pettis & Shimanuki, 2000), but a minimum depth rather than a preferred depth is likely (Meikle & Diaz, 2012). In our study we assumed a depth of 10 cm to estimate pupal performance. In (semi‐)desert top soils moisture levels could fluctuate and even dry out completely. This could affect the chances to complete the life cycle and thus establishment of populations of SHBs under these extreme conditions.

Small hive beetles can survive temporarily unsuitable environmental conditions due to the thermoregulatory capacity of honey bee colonies (Schäfer, Ritter, Pettis, & Neumann, 2011). SHBs have been shown to survive in honey bee winter clusters for several months. However, SHBs cannot maintain populations close to the temperate climatic limits of beekeeping (Neumann et al., 2016) as we found these areas to be unsuitable for pupation. Therefore, within the predicted current and future range expansion of SHB presented here, honey bee colonies will inevitably be present and thus potentially exposed to this pest species. Moreover, SHBs infest colonies of other social bees as well (Cervancia et al., 2016; Greco et al., 2010; Spiewok & Neumann, 2006) and may also use solitary bees as hosts (Gonthier et al., 2019). Even though the impact and its magnitude are unknown for many social bees, the predicted invasion risk could have a broad impact on this important group of pollinators and the ecosystem service they provide. While it is well known that European honey bee subspecies are more susceptible to SHB infestations compared to African ones, probably due to quantitative differences in a range of defence behaviours (Ellis & Hepburn, 2006; Neumann & Elzen, 2004), there are no data available for other honey bee or other bee species. This notion vastly expands the impact of SHB as an invasive species, having potential consequences for species already declining, such as certain bumble bees (Meeus, Brown, Graaf, & Smagghe, 2011; Potts, Roberts, et al., 2010).

With our mechanistic niche models, we assessed critical dimensions of the fundamental niche of SHBs, which allowed a quantification of the direct impacts of ambient climate on pupal performance outside the host colony. This provides a first basis for risk assessment at the global scale, but local realisations might still be modified to some extent by factors impacting other life stages, e.g. by dispersal limitation, different factors affecting establishment and epidemiology or local soil physics. The first step of successful invasion is the introduction (Richardson et al., 2000), which can be by natural spread of the beetle alone or via the bees, once it arrives in a new area, or by anthropogenic movement of colonies or bee products (Lounsberry et al., 2010; Neumann et al., 2016). Given the well‐documented role of global trade and movement of bees and bee products for the spread of SHBs (Neumann et al., 2016; Ouessou Idrissou et al., 2019), introduction seems not to be a limiting factor for the invasion success of SHB.

For successful establishment and further spread, the presence and density of known and new hosts as primary and alternative food sources can be important (Gonthier et al., 2019; Schäfer et al., 2019). In the native range of sub‐Saharan Africa, feral honey bee populations show higher densities compared to other regions in the world (Jaffé et al., 2010). However, in Asia, Europe and the United States managed honey bee populations are much more prolific and colonies are concentrated in apiaries, probably affecting invasion dynamics of SHBs accordingly (Neumann et al., 2010; Spiewok, Duncan, Spooner‐Hart, Pettis, & Neumann, 2008). Especially when considering the good flight ability of adult SHBs (Neumann, Hoffmann, Duncan, Spooner‐Hart, & Pettis, 2012), such high host densities likely benefit the establishment and further spread of SHB in addition to movement of colonies or bee products (reviewed by Neumann et al., 2016), in particular in the absence of early detection systems (Schäfer et al., 2019). Furthermore, the host could also be affected by climatic conditions, indirectly altering the conditions for SHB in different life stages than the pupa. The effect of climate change on honey bees, however has not been well studied, but is likely to include effects on phenology and survival (Le Conte & Navajas, 2008). For instance, elevated temperatures can expedite the onset of brood rearing in honey bee colonies in temperate climate zones (Nürnberger, Härtel, & Steffan‐Dewenter, 2018). Furthermore, the availability of food sources and the opportunity to forage could change (Le Conte & Navajas, 2008), thereby affecting brood rearing cycles. SHB follows honey bee phenology (Lundie, 1940) and likely seizes to reproduce in winter, as only adults are found in winter clusters (Schäffer et al., 2011). The availability of brood in winter could trigger the onset of reproduction in SHB, but seasonality of SHB reproduction in temperate climatic zones and the relation to its host has not yet been studied (Neumann et al., 2016). While our study provides a first identification of suitable areas and potential severity of SHB under current and future conditions, which can be used for precautionary management plans, local realisations, e.g. by identifying the relative importance of soil conditions, dispersal and anthropogenic translocations or indirect impacts of global change on the host species, still warrant further investigations.

Generally, our study further highlights the urgent need for slowing down the global spread of SHBs (Schäfer et al., 2019), until better mitigation options are available. The results provide a science‐based approach in support of strategic management of this invasive species as measures can be taken where they are deemed fit for now or for the future. Past invasions of SHB have resulted in an economic deficit to the beekeeping industry (Neumann & Elzen, 2004; Rhodes & McCorkell, 2007). The focus of management strategies should therefore firmly be on detection in the early stages of invasion (Hulme, 2009) by focussing on the global trade in bees and bee products (Neumann et al., 2016; Ouessou Idrissou et al., 2019). Options for which are provided by improving the international trade regulations (Lecocq et al., 2016).

In conclusion, our study shows for the first time an increased global invasion risk of a honey bee pest due to global warming. As managed honey bee populations and many wild bee species are either suffering from unsustainably high colony losses (Neumann & Carreck, 2010) or decline (Potts, Roberts, et al., 2010), it appears prudent to further investigate the interplay between climate change and biological invasions in the context of bee health. Furthermore, our study merits further investigation of the potential risks posed by other alien, invasive honey bee pests like Varroa destructor and Tropilaelaps sp. under climate change scenarios. But even native pests and pathogens should not be excluded from such analysis (Le Conte & Navajas, 2008). In more general terms, we need a better understanding of the impact of climate change on biological invasions and the impacts on ecosystem services (Knight et al., 2018; Phillips et al., 2018).

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTION

B.C., O.S. and P.N. developed the idea and wrote the manuscript. B.C. compiled peer‐reviewed data, collected presence data and contributed additional laboratory experiments. O.S. designed and built models, performed model analyses and analysed data.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.