Predicting Potential Global Distribution of the Invasive Species Aethina tumida Murray (Coleoptera: Nitidulidae), and Its Natural Enemies Steinernema carpocapsae (Weiser, 1955)

1. Introduction

Global climate change, particularly rising temperature and shifting precipitation patterns, is profoundly reshaping the geographical distribution and suitable habitats of species, driving range shifts toward higher latitudes and elevations [1,2]. Temperature serves as a critical environmental determinant regulating insect distribution, directly constraining survival, developmental rates, and reproductive potential. Within optimal thermal thresholds, warming can enhance insect growth and population proliferation; however, once tolerance limits are exceeded, heat stress induces physiological disorders, reduces reproductive fitness, and can precipitate population collapse [3]. Concurrently, invasive insects substantially alter the population dynamics and spatial distribution of indigenous species through competitive exclusion, predation, or pathogen transmission, thereby disrupting established interspecific interaction networks and precipitating irreversible degradation of ecosystem structure and function. Of particular concern is the marked synergism between climate change and anthropogenic activities—including international trade, land-use modification, and transportation networks. Elevated temperatures expand the potential habitat range of invasive species, while globalization has markedly increased the probability of incidental insect introduction and establishment via commodity movement and transport infrastructure, consequently accelerating their spread and exacerbating biodiversity loss and ecosystem disequilibrium [4,5]. Invasive alien species frequently modify the population size and distribution patterns of native species, disrupting regional species interrelationships and subsequently inducing irreversible alterations in the structure and function of invaded ecosystems. Such incursions not only cause inflict substantial economic damage upon agriculture, forestry, animal husbandry, and fisheries, but also compromise environment safety and public health [6]. Consequently, comprehensive analysis of the complex interactions between invasive insects and local biotic communities as well as abiotic environmental factors is imperative to elucidate the ecological drivers underlying outbreak dynamics, and to develop precision monitoring and early-warning systems alongside environmentally sustainable management technologies, thereby systematically mitigating the ecological and economic risks posed by invasive insects [7].

The small hive beetle, Aethina tumida Murray (Coleoptera: Nitidulidae), is a destructive invasive pest of honey bees originally distributed in sub-Saharan Africa, where it is widespread in tropical and subtropical regions [8,9]. It has been listed as one of the six important bee pathogens by the World Organization for Animal Health (WOAH) (https://www.woah.org/en/what-we-do/animal-health-and-welfare/animal-diseases/) [10,11,12]. Propelled by the rapid expansion of the beekeeping industry and increased trade in bee products, this pest has disseminated to North America, South America, Asia, and Oceania, with reports now documented in nearly twenty countries worldwide [13,14,15]. Both adult and larval stages feed upon bee brood, honey, and pollen, excavating comb structures and compromising hive integrity. Larval excrement induces honey fermentation, generating foul odors that precipitate colony absconding [16,17]. Infestations result in discolored，fermented honey emitting a distinctive odor resembling of rotten oranges. When nest structures and hive covers sustain damaged and fermentation, honey may effervesce bubble and leak from the colony. Larvae frequently deposit viscous, malodorous residues that induce bee abandonment. Adult beetles exhibit predilection for bee eggs and larvae, severely compromising colony development and potentially causing colony collapse, absconding behavior, or mortality [18,19]. Given the photophobic behavior of adults, they seek refuge in corners and crevices upon hive inspection. Detection of adults at the hive bottom necessitates vigilance regarding potential larval damage. Early detection proves challenging because larvae remain concealed within sealed cells [20]. A. tumida consumes all hive products and disperses through multiple mechanisms, including autonomous flight, transportation of infected colonies, and human-mediated assistance. Primary dispersal pathways encompass migratory beekeeping operations, include migrating bee colonies, beehives and beeswax, soil attached to various used for import-export packaging, and parasitism in circulating fruits and vegetables [21,22,23,24]. Furthermore, bumblebees have been identified as susceptible hosts, and beyond honey bees, A. tumida can parasitize alternative bee populations including stingless bees (Apidae: Meliponini) [25,26], bumblebees (Bombus spp.) [27], and solitary bees [28]. Additionally, A. tumida functions act as a vector for transmission of honey bee pathogens, including Pasteurella larvae [29,30] and viruses such as deformed wing virus and sac brood virus [31]. Stenernema carpocapsae (Weiser, 1955), an entomopathogenic nematode characterized by broad insecticidal spectrum, facile cultivation, low resistance potential, and safety to humans, livestock, and the environment, has emerged as one of the most promising biological control agents [32,33].

Species distribution models (SDMs) constitute effective tools for predicting the spatial distribution of species’ suitable habitats. Commonly used SDMs include BIOCLIM, GARP, MaxEnt, CLIMEX, and DOMAIN, which have been extensively applied across diverse research domains [34,35]. MaxEnt model currently ranks among the most widely utilized species distribution modeling approaches. This model predicts potential species distribution range based on documented occurrence records and environmental variables [36,37]. CLIMEX represents a semi-mechanical modeling methodology emphasizing the ecological and physiological responses of species to their environmental niche, whereas MaxEnt emphasizes the statistical relationship between species occurrence localities and environmental variables. Drawing upon documented geographical distributions and biological characteristics data, CLIMEX simulates potential suitable habitat ranges under climatic conditions, that optimally correspond to actual colonization patterns within documented distribution ranges, thereby predicting potential suitable habitats within target regions [38].

This study employed MaxEnt and CLIMEX models to predict the global potential distribution for the invasive species A. tumida and its natural enemy S. carpocapsae under current and future climate conditions. The primary objectives were to: (1) identify the key environmental factors influencing the potential geographical distribution of both species; (2) forecast potential distribution patterns and dynamic range shifts under climate change scenarios; and (3) establish a theoretical framework to inform prevention, monitoring, and sustainable management strategies.

2. Materials and Methods

2.1. Occurrence Data

Occurrence data were obtained from the following sources: for A. tumida, the Global Biodiversity Information Facility (GBIF, http://www.gbif.org/, https://doi.org/10.15468/dl.bajcbm) and relevant published research literature [12]. For S. carpocapsae, the Global Biodiversity Information Facility (GBIF, http://www.gbif.org/, https://doi.org/10.15468/dl.8nqwrb), the Centre for Agriculture and Bioscience International (CABI, https://plantwiseplusknowledgebank.org/doi/full/10.1079/pwkb.species.51706), and previously published literature [39,40,41]. Data preprocessing was conducted using the dismo package in R 4.5.2 was to eliminate duplicate records, geospatially erroneous points, and observations lacking geographic information. To mitigate sampling bias and spatial autocorrelation, the Spatial Rarefy Occurrence Tool in SDM Toolbox 2.5 was employed to thin the dataset, applying a minimum distance threshold of 5 km between retained points [42]. Finaly, the final modeling dataset comprised 225 for A. tumida and 94 for S. carpocapsae, respectively. Additionally, host distribution data were retrieved from the Global Biodiversity Information Facility Database (GBIF, https://www.gbif.org, https://doi.org/10.15468/dl.4hsz5a), providing critical context for analyzing the relationship between host availability and target species distribution (Figure 1).

2.2. Bioclimate Data

Current bioclimatic variables were obtained from the World Climate Database (http://www.worldclim.org/) at a spatial resolution of 2.5 arc-minutes. This dataset comprises 19 bioclimatic variables (bio1-bio19), representing minimum, maximum, and mean values of monthly, seasonal, and annual ambient temperatures for the baseline period 1970–2000 [43] (Appendix A, Table A1). Climate data for CLIMEX modeling were obtained from the CliMond dataset (https://www.climond.org/) with a spatial resolution of 30 arc-minutes based on the Special Report on Emission Scenarios (SRES). This dataset encompasses long-term monthly averages of minimum and maximum temperatures, together with monthly mean relative humidity at 09:00 (RH 0900) and 15:00 (RH 1500), spanning the period1961-1990and centered on 1975 [44].

Future bioclimatic variables were derived from the Coupled Model Intercomparison Predict Phase 6 (CMIP6) models [45]. Three Shared Socioeconomic Pathways (SSPs)—SSP126, SSP245, and SSP585—were selected to represent low, medium, and high greenhouse gas emission trajectories, respectively [46], across four future periods: the 2030s (2021–2040), 2050s (2041–2060), 2070s (2061–2080), and 2090s (2081–2100). For CLIMEX model, the A1B and A2 scenarios were employed for the 2030s, 2050s, 2070s, 2090s, and 2100s [47]. The A1B scenario characterizes rapid global economy integration, diversified and balanced energy development, and the progressive substitution of fossil fuels by non-fossil energy sources as the dominant primary energy component, predicting a global mean surface temperature increase of 1.4–3.8 °C by 2100 (best estimate: 2.4 °C). The A2 scenario depicts a regionally heterogeneous development pattern characterized by sustained global population growth and fragmented, decentralized economic expansion, assuming atmospheric CO₂ concentrations of 846 ppm and predicting approximately 6 °C warming by the end of this century [48].

Environmental variables exhibit substantial spatial correlation, which may compromise predictive accuracy and induce model overfitting [49]. To address multicollinearity, pearson correlation analysis was performed using IBM SPSS Statistics 25 (Appendix B, Figure A1), and variables exhibiting the strongest autocorrelation (|r| > 0.8) were excluded [50,51]. The final variable set comprised five bioclimatic predictors for A. tumida (bio2, bio10, bio11, bio12, bio18) and five for S. carpocapsae (bio2, bio11, bio17, bio18, bio19).

2.3. MaxEnt Model

The MaxEnt model (3.4.4; https://biodiversityinformatics.amnh.org/open_source/maxent/) implements maximum entropy theory to estimate species’ potential distributions by contrasting occurrence localities against randomly sampled background points. Model parameterization followed established protocols: 75% of occurrence records were allocated selected for model training, with the remaining 25% reserved for independent validation; regularization multiplier was set to 1; repetitions was fixed at 10 runs; maximum iterations were constrained to 500; and 10000 random background points were generated. Variable importance was assessed via the Jackknife method, and response curves were generated to characterize species–environment relationships. Model output exported as continuous raster data in logistic format, with all remaining parameters retained at default values [52,53]. The result was categorized into four classes using the reclassification tool in ArcGIS 10.8: unsuitable (0–0.2), marginally suitable (0.2–0.4), moderately suitable (0.4–0.6), and highly suitable (0.6–1) [54]. To analyze the spatiotemporal dynamics of species distribution, continuous outputs were converted to binary presence/absence maps using the thresholding tool in the SDM Toolbox, from which range shifts were classified as range expansion, no change, or range contraction [55,56].

Model performance was evaluated using the area under the receiver operating characteristic curve (AUC) and true skill statistics (TSS) [52]. AUC ranges from 0 to 1 and integrates sensitivity and specificity; values >0.9 indicate excellent performance, 0.8-0.9 good, and <0.7 poor [57,58]. TSS is calculated using the validation dataset and is not affected by the size of the validation dataset, ranges from -1 to 1, with values approaching 1 indicating optimal agreement between predictions and observations; TSS >0.75 denotes excellent model accuracy [59].

Predictive accuracy in MaxEnt is primarily governed by two important parameters: the regularization multiplier (RM) and feature combination (FC). To optimize model these parameters and mitigate overfitting, the R package ENMeval was employed [60]. Feature classes evaluated comprised linear (L), product (P), quadratic (Q), threshold (T), and hinge (H) functions [61]. Eight FC combinations were assessed as candidate models: L, LQ, LQH, LQHP, LQHPT, LQP, QHP, and QHPT. RM was incrementally varied from 0.5 to 4.0 in 0.5-unit steps [62]. Model selection was based on the Akaike information criterion corrected for small sample sizes (AICc), with the model exhibiting the lowest delta AICc (ΔAICc=0) designated as optimal [63,64]. Following optimization, the selected parameters were FC=LQHPT and RM=1 for A. tumida, and FC=LQH and RM=1 for S. carpocapsae, both with ΔAICc=0 (Appendix B, Figure A2).

2.4. CLIMEX Model

CLIMEX v.4.0 (Hearne Scientific Software, Australia) was implemented using the “Compare Locations” function to construct simulation models [65,66]. Key biological parameters incorporated included: developmental temperature thresholds (DV0, DV1, DV2, DV3); soil moisture thresholds (SM0, SM1, SM2, SM3); degree-days per generation (PDD); cold stress parameters (threshold temperature, TTCS; accumulation rate, THCS); heat stress parameters (threshold temperature, TTHS; accumulation rate, THCS); dry stress parameters (threshold temperature, SMDS; accumulation rate, HDS); and wet stress parameters (threshold soil moisture, SMWS; accumulation rate, HWS) [67]. These parameters were calibrated iteratively against known distribution records, with unlisted parameters retaining system default values (Appendix A, Table A2). Through repeated model runs, parameter refinement, and cross-validation against occurrence data, the definitive parameter set was established. Calibration was validated by ensuring modeled distributions encompassed all documented occurrence localities, with final outputs mapped in ArcGIS to visualize suitable ranges [68]. Model outputs were expressed as Ecoclimatic Index (EI) values, ranging from 0 to 100, where EI = 0 indicates the region is unsuitable for species survival and EI approaching 100 indicates that the climate and environment in the region are more suitable [69,70]. The value was classified as: unsuitable (EI=0), marginally suitable (0<EI<10), moderately suitable (10<EI<20), and highly suitable (EI>20) [71].

3. Results

3.1. Model Accuracy Evaluation and Host Availability

Following 10 replicate runs, the MaxEnt models demonstrated robust predictive performance, with A. tumida achieving AUC=0.959 and TSS=0.815, and S. carpocapsae AUC=0.900 and TSS=0.842. These indicates excellent model accuracy and high reliability for predicting potential suitable habitats (Appendix B, Figure A3).

Percentage contribution analysis revealed that annual precipitation (bio12, 43.1%), mean temperature of warmest quarter (bio10, 24.5%), and mean temperature of coldest quarter (bio11, 23.6%) were the predominant determinants of A. tumida distribution, followed by collectively accounting for 91.2% of explained variance. For S. carpocapsae, the most important environmental factors were mean temperature of coldest quarter (bio11, 58.3%), precipitation of driest quarter (bio17, 22.2%), and precipitation of coldest quarter (bio19, 13.9%) emerged as the most influential variables, contributing 94.4% cumulatively (Appendix B, Figure A4). Jackknife regularization tests corroborated these findings, identifying bio12 and bio11 as the variables with highest predictive power for A. tumida and S. carpocapsae, respectively, both exhibiting training gains exceeding 0.82 (Appendix B, Figure A5). Response curve illustrated species–environmental variables with high clarity: optimal conditions A. tumida occurred at bio2 = 8-15 °C, bio10 = 18-30 °C, bio11 = -10 to 25 °C, bio12 peaking at 1010mm (suitability = 0.7), and bio18 = 0-800mm (Appendix B, Figure A6). For S. carpocapsae, suitable ranges encompassed bio2 = 8-16 °C, bio11 = -10 to 20 °C, bio17 = 0–1000mm, bio18 = 0-800mm, bio19 = 0-600mm (Appendix B, Figure A7).

3.2. Potential Distribution of A. tumida and Host Availability

It is worth noting that the potential distribution of A. tumida is basically overlaps with the spatial distribution of host plants. The results indicate that host plants are widely distributed across all continents (Figure 1 and Figure 2).

3.3. Potential Distributions Under Climate Conditions Using MaxEnt

The MaxEnt model is used to predict potential shifts in the global geographic distributions of A. tumida and S. carpocapsae under current and future SSP scenarios (SSP126, SSP245, SSP585) across four different periods (2030s, 2050s, 2070s, and 2090s) (Figure 2, Figure 3, Figure 4 and Figure 5).

Under current climate conditions, highly suitable habitat for A. tumida occupied 3.75 million km2 (13.03% of the total suitable area), concentrated predominantly in the United States. Moderately suitable habitat encompassed 8.39 million km2 (29.15%), distributed across eastern Argentina, Southern Brazil, Angola, Zambia, and parts of southern Africa, whereas marginally suitable habitat comprised 16.64 million km2 (57.82%), distributed in central region of South America, southern Africa, and southern Asia. Under future climate scenarios, divergent distributional trajectories emerge: SSP245 (2050s) and SSP585 (2070s) have marked expansion into North America, Europe and central Australia, whereas SSP585 (2050s) indicates pronounced contractions across Asia, Africa, and South America.

For S. carpocapsae, highly suitable habitat totaled 10.09 million km2 (14.96% of total suitable area), with distributions in southern China, the United States and parts of Europe. Moderately suitable habitat covered 21.92 million km2 (32.49%), concentrated in the United States, portions of South America, Europe, and Asia, whereas marginally suitable habitat extended across 35.46 million km2 (52.56%), encompassing all continents. Future projections reveal complex spatiotemporal dynamics: substantial range expansion into Asia, North America, and Africa is anticipated under SSP126 (2090s), SSP245 (2030s), and SSP585 (2070s), yet significant contractions across these same regions are projected under SSP585 by the 2030s.

3.4. Potential Distributions Under Climate Conditions Using CLIMEX

CLIMEX employed to predict the potential global distribution shifts for A. tumida and S. carpocapsae, under current and two climate scenarios (A1B, A2) across five temporal horizons (2030s, 2050s, 2070s, 2090s and 2100s) (Figure 2, Figure 3, Figure 6 and Figure 7).

Under current climate conditions, the potential suitable habitat of A. tumida encompassed South America, southern Asia, and the majority of the Africa. Under A1B and A2 scenarios, total suitable habitat exhibited incremental expansion of 0.94% and 2.38%, respectively, by 2100s. However, this aggregate stability masked substantial internal reconfiguration: highly suitable habitat contracted markedly by 22.73% (A1B) and 33.08% (A2), whereas moderately suitable habitat expanded by 78.65%(A1B) and 131.96% (A2). Similarly, marginally suitable habitat is projected to increase by 45.78% (A1B) and 62.76% (A2), respectively, by 2100.

Under current climatic conditions, suitable habitat of S. carpocapsae is distributed across all continents. Climate change predictions indicate substantial restructuring of this distributional pattern. Under both A1B and A2 climate scenarios, highly suitable habitats are predicted to contract markedly, whereas moderately and marginally suitable habitats are expected to expansion. Specifically, by the 2100s, highly suitable habitat is anticipated to decrease by 23.42% (A1B) and 22.54% (A2). In contrast, moderately suitable habitat exhibits pronounced expansion by 91.74% under A1B (2100s) and 86.03% under A2 (2090s). Similarly, marginally suitable habitat is projected to increase by 10.29% (A1B) and 21.27% (A2), respectively, by 2100.

3.5. Combined Prediction Maps of the Two Models of A. tumida and S. carpocapsae Under Current Climate Conditions

To enhance predictive robustness, outputs from MaxEnt and CLIMEX were integrated under current climatic conditions, with spatial overlap between model predictions interpreted as high-confidence suitable habitat. Both models concurred in predicting pantropical and temperate distribution potential across all continents. Notably, CLIMEX predictions encompassed the entire spatial domain predicted by MaxEnt (Figure 2).

4. Discussion

This study employed MaxEnt and CLIMEX models to predict potential suitable habitats for A. tumida and S. carpocapsae under current and future climate change. Following parameter optimization (FC=LQHPT, RM=1 for A. tumida; FC=LQH, RM=1 for S. carpocapsae), both models achieved AUC and TSS values exceeding 0.8, indicating predictive accuracy and high reliability [72,73]. The relative importance of bioclimatic variables differed markedly between species. A. tumida distribution was predominantly constrained by annual precipitation (bio12), mean temperature of warmest quarter (bio10) and mean temperature of coldest quarter (bio11), corroborating previous findings [12]. Conversely, S. carpocapsae suitable habitat was principally governed by mean temperature of coldest quarter (bio11) and precipitation of driest quarter (bio17). Temperature and precipitation thus emerge as critical climatic limiting factors for both species, consistent with established ecological theory [74]. Temperature exerts direct physiological control over A. tumida developmental rates and survival across life stages. The lower thermal thresholds for egg, larval, and pupal development range from 10.0–13.5 °C [75], whereas temperatures exceeding 35℃ severely compromise egg viability. Under favorable thermal regimes, A. tumida can complete up to six generations annually. Pupae survival and developmental duration are temperature-dependent, with optimal conditions prolonging development yet enhancing survival [76]. Notably, A. tumida remarkable cold tolerance, capable of overwintering within hives at -40℃, though it cannot persist outside hive structures in northern latitudes, restricting reproduction to summer months. Humidity further modulates population dynamics: relative humidity below 50% facilitates egg hatching, and larval-to-pupal development typically requires 10-14 days, with pupation rate of 92%-98%recorded in moist soil [21]. As a thermophilic and hygrophilic pest, A. tumida is poised to expand its distributional range under global warming scenarios, posing persistent threats to apicultural systems worldwide threat to apicultural system worldwide [77].

Under current climate conditions, both models indicate that suitable habitat for A. tumida was concentrated in South America, southern Africa, and southern Asia, whereas S. carpocapsae exhibited a broader, near-cosmopolitan potential distribution. Comparative analysis revealed that CLIMEX predicts encompassed more extensive suitable ranges than MaxEnt for both species, with CLIMEX outputs effectively subsuming MaxEnt predictions. These discrepancies reflect fundamental methodological differences between modeling approaches. MaxEnt relies primarily on species occurrence records and environmental covariates, offering flexibility to incorporate diverse predictor variables including edaphic and topographic factors. CLIMEX, by contrast, integrates species distribution data with eco-physiological parameters, explicitly modeling four stress indices (hot, cold, dry, and wet) to characterize species-climate relationships. Ensemble integration of both model outputs consequently enhances predictive confidence and reduces algorithm-specific uncertainty [68].

A. tumida is a nest pest of the western honey bee, Apis mellifera L. [78], which a species indigenous to Europe, Africa, and the Middle East that has achieved cosmopolitan distribution through anthropogenic dispersal. Intensified international trade and bee stock exchange have facilitated A. tumida range expansion, inflicting severe damage on apicultural operations [14]. Beyond structural hive degradation, infestations precipitate colony depopulation and substantial economic losses. S. carpocapsae, as a key natural enemy of A. tumida, demonstrates considerable potential for biological control application, effectively suppressing pest population and furnishing critical technical support for sustainable, environmentally benign management strategies.

This study incorporated 19 bioclimatic variables to predict suitable habitats for A. tumida and its natural enemy S. carpocapsae. Nevertheless, edaphic factors such as soil temperature and moisture, which may significantly influence on species distribution [79]. Furthermore, the models did not account for critical biotic and abiotic factors—interspecific competition, anthropogenic disturbance, and dispersal barriers that constrain the ability to explain actual invasion dynamics. The beehive microclimates and other environmental features from macroclimatic models necessarily reduced predictive precision.

Future research should integrate species-specific biological characteristics, beehives microclimate, and wild habitats characteristics, incorporating environmental covariates such as soil physicochemical properties and topographic factors, while supplementing biotic and abiotic factors including interspecific competition and human activities. Such comprehensive model optimization would provide robust theoretical foundations for biosecurity management.

Global agricultural systems exhibit substantial dependent on bee-mediated pollination services, particularly for oilseeds crops, forage plants, and fruits and vegetables production, conferring significant economic and ecological value. Invasive pests disseminate across international borders via trade and transportation networks, posing severe threat to apicultural industries and ecological security. Port quarantine protocols should prioritize high-risk vectors including bee products, wooden packaging materials, and soil-adhered horticultural commodities, with concomitant strengthening of access management and early-warning systems. For regions where establishment has already occurred, adaptive zonal establishment should be implemented: enhanced monitoring and emergency response capacity in highly suitable habitats; establishment of buffer zones in peripheral habitats; and development of transregional joint prevention and control mechanisms.

5. Conclusions

The potential suitable distribution for A. tumida predominantly concentrated in South America, Africa, and Asia. Under scenarios of continued global environmental changes and intensified anthropogenic activities, this species exhibits a pronounced poleward range expansion trajectory. Given its substantial invasive potential and associated ecological risks, urgent implementation of enhanced prevention and control measures is imperative to mitigate threats to ecosystems and agricultural production systems. Concurrently, S. carpocapsae has achieved near-cosmopolitan distribution owing to its broad environmental tolerance and effective dispersal mechanisms, underscoring the necessity for sustained monitoring and targeted management strategies for this biological control agent.