Data and code from: Mapping the global range of stingless bees: Insights into their ecological niches and environmental tolerances
Data files
Sep 21, 2026 version files 305.45 MB
-
Data.zip
47.61 MB
-
Outputs.zip
257.82 MB
-
README.md
10.11 KB
-
Scripts.zip
19.05 KB
Abstract
Aim
Stingless bees are the most species-rich group of eusocial bees, predominantly found in tropical and subtropical regions of the world. However, despite their diversity and importance, significant gaps remain in our understanding of their distribution patterns and the factors driving them, as these aspects have yet to be thoroughly investigated on a global scale. In this study, we aimed to fill this gap by mapping global habitat suitability for stingless bees.
Location
Global
Taxon
Stingless bees (Apidae:Apinae:Meliponini)
Method
We compiled over half a million spatiotemporal occurrence records from 361 species of stingless bees at an unprecedented scale from both public and private databases optimised using a recently published data-cleaning workflow. Through a species distribution model approach, we assessed the relationship between stingless bee distribution, diversity, and environmental factors, ultimately characterising their ecological niche and producing the first habitat suitability maps for the bee tribe Meliponini.
Result
Our results reveal significant discrepancies in occurrence records across the Neotropical, the Afrotropical, and the Indo-Malayan, Papuasian, and Australian regions, as well as among taxa. Habitat suitability hotspots were identified in tropical forests and were associated with relatively high temperatures and humidity throughout the year. Projection maps further suggested that some temperate, xeric, and savannah regions could also be suitable habitats for stingless bees. Our results rescale the known environmental niche of stingless bees, demonstrating that suitable habitats exist well beyond the tropical and subtropical boundaries typically described in the literature.
Main conclusion
This study presents a novel, model-based assessment of stingless bee distribution, revealing that it extends beyond tropical forests into drier, colder, and temperate regions more broadly than previously recognized. These findings highlight substantial gaps in ecological and distributional data, emphasizing the need for taxonomic revision, targeted fieldwork, digitisation of collections, and foundational research on stingless bee ecology and evolution. Addressing these gaps will improve predictions of species responses to environmental change, inform conservation strategies, and advance understanding of stingless bees as essential providers of ecosystem services.
Dataset DOI: 10.5061/dryad.kkwh70sfw
This README outlines the data files, scripts, and figures used in the analysis, allowing for replication of results and further exploration. Through a species distribution model approach, we assessed the relationship between stingless bee distribution, diversity, and environmental factors, ultimately characterising their ecological niche and producing the first habitat suitability maps for the bee tribe Meliponini.
The repository is organized into the following directories (zip files):
Data: Contains the public data used in the analysis
Scripts: Contains R scripts necessary for reproducing the analysis.
Outputs: Contains output of models for each species
1. Data Files
1.1 Data.zip
UnCleanStingless_2023-01-25
Contains the occurrence records extracted from public repositories using the beeDC workflow without any filtering of the flagged occurrences on September 25th 2023.
cleanStingless_2023-09-25.csv
Contains the occurrence records extracted from public repositories using the beeDC workflow filtering the flagged occurrences on September 25th 2023.
Both Unclean and Clean datasets contain geographic occurrence records of stingless bees (Apidae: Apinae: Meliponini). Each row represents an occurrence record and contains taxonomic identification, geographic coordinates, and temporal information where available.
The data are provided as a semicolon-delimited (;) text/CSV file.
Data dictionary
| Variable | Description | Units/Values |
|---|---|---|
| ID | Unique occurrence record identifier. | Text |
| BinomialName | Full scientific species name. | Text |
| Family | Taxonomic family | Apidae |
| Subfamily | Taxonomic subfamily | Apinae |
| Tribe | Taxonomic tribe | Meliponini |
| Genus | Taxonomic genus. | Text |
| Species | Specific epithet. | Text |
| Lat | Latitude of occurrence. | Decimal degrees (°); −90 to +90 |
| Long | Longitude of occurrence. | Decimal degrees (°); −180 to +180 |
| Year | Year of occurrence/collection. | YYYY |
| Month | Month of occurrence/collection. | 1–12 |
| Day | Day of occurrence/collection. | 1–31 |
| Date | Complete occurrence/collection date. | ISO 8601; YYYY-MM-DDTHH:MM:SSZ |
sp_occ_grid.csv
Contains expert-derived species occurrence records represented as grid-cell centroids. In the Species R script, the original expert dataset (data_expert.csv) is converted to grid-cell centroids for computational analyses. However, some data contributors do not permit the sharing of exact geographic coordinates. Therefore, for reproducibility, sp_occ_grid.csv should be used as the input file instead of data_expert.csv.
Data dictionary
| Variable | Description | Units/Values |
|---|---|---|
| X | Longitude of occurrence. | Decimal degrees (°); −180 to +180 |
| Y | Latitude of occurrence. | Decimal degrees (°); −90 to +90 |
| Species | Full scientific species name. | Text |
Grid/
This folder contains the 25 km² grid shapefile used for spatial analyses. The shapefile consists of the following four files:
ShpGrid25km_world_Clip.shp = geometry of the grid cells
ShpGrid25km_world_Clip.shx = spatial index
ShpGrid25km_world_Clip.dbf = attribute data
ShpGrid25km_world_Clip.prj = coordinate reference system information
Environmental predictors.zip
This compressed folder contains the 35 raster (.tif) files representing the environmental predictors used in this study. The rasters were cropped to exclude Antarctica and standardized (scaled) prior to analysis.
1.2 Scripts.zip
This folder contains the scripts necessary to the analysis described in the manuscript:
Species.R
Main script that is divided into 9 sections: 1) list of packages and custom function 2) Data import & Cleaning/filtering occurrences 3) create datasets 4) Descriptive plots 5) Create/import grid 6) Import predictors 7) Geographical sampling bias 8) SDMs - Loop for each species 8) Plot prediction maps & predictors importance.
create_grid.R
Custom function used in Species.R to create a global 25km² grid.
1.3 Outputs.zip
This folder contains the outputs produced by each model using the Species.R script. For each species, several files are generated following the same naming structure:
Genus_Species_File_Name_RunNumber.extension
RunNumber corresponds to the repetitions of the same model computation, which are performed to account for uncertainty.
Note that figures in the manuscript were produced based on a subset of species based on the model performance, see Script.R and the method section of the manuscript.
Null models
Genus_Species_mod_perform_enmeval_null_runXX.csv
MaxEnt null model evaluation metrics. Columns:
| Column | Description | Unit / value |
|---|---|---|
Codes |
Unique model identifier | Character; equal to fc.rm |
fc |
Feature classes used in MaxEnt | Character; L or LQ |
rm |
Regularization multiplier | Numeric; 0.5-5 |
tune.args |
MaxEnt tuning specification | Character |
auc.train |
Training AUC | 0–1 |
cbi.train |
Training Continuous Boyce Index | −1 to 1 |
auc.diff.avg |
Mean difference between training and validation AUC | AUC |
auc.diff.sd |
SD of training–validation AUC difference | 0-1 |
auc.val.avg |
Mean validation AUC across replicates | 0–1 |
auc.val.sd |
SD of validation AUC across replicates | 0-1 |
cbi.val.avg |
Mean validation Continuous Boyce Index | −1 to 1 |
cbi.val.sd |
SD of validation CBI across replicates | 0-1 |
or.10p.avg |
Mean omission rate using the 10th-percentile training threshold | 0–1 |
or.10p.sd |
SD of the 10th-percentile omission rate | 0-1 |
or.mtp.avg |
Mean omission rate using the minimum training presence threshold | 0–1 |
or.mtp.sd |
SD of the minimum training presence omission rate | 0-1 |
AICc |
Corrected Akaike Information Criterion | AICc score |
delta.AICc |
Difference between model AICc and minimum AICc | ΔAICc |
w.AIC |
AICc model weight | 0–1 |
ncoef |
Number of model coefficients | Count |
Species |
Species modeled | Character |
or.10p.avg.round |
Rounded mean 10th-percentile omission rate | 0–1 |
auc.val.avg.round |
Rounded mean validation AUC | 0–1 |
Null.AUC.95 |
95th percentile of the null AUC distribution | 0–1 / NA |
Genus_species_null_model_runs_runXX.pdf
Histogram representing the reference distribution of validation AUC values. Vertical lines indicate the 1st, 5th, 50th, 95th, and 99th percentiles, while the red line indicates the empirical AUC value.
Model results
Genus_species_mod_perform_enmeval_runXX.csv
MaxEnt null model evaluation metrics with the same variable than for null models.
Genus_species_var_importance_enmeval_run10.csv
Relative contribution (reflecting the influence during model training) and permutation importance (measuring prediction changes with randomised variable values) for each predictor included in the species-specific model.
Genus_species_response_curve.pdf
Line plots of the response curves for each predictor included in the species-specific model.
Model Predictions
Genus_species_prediction_present_wide_runXX.pdf /.tif
Habitat suitability projected at the global scale.
Genus_species_median_prediction_present_wide_runXX.tif
Median habitat suitability projected at the global scale. The value represents the median across all 10 model iterations.
Genus_species_sd_prediction_present_wide_runXX.tif
Standard deviation of habitat suitability projected at the global scale. The value in each grid cell represents the standard deviation across the 10 model iterations.