Research Article
Research Article
Niche modeling of bumble bee species (Hymenoptera, Apidae, Bombus) in Colombia reveals highly fragmented potential distribution for some species
expand article infoLaura Rojas-Arias, Daniel Gómez-Morales§, Stephanie Stiegel|, Rodulfo Ospina-Torres§
‡ National Autonomous University of Mexico, Mexico City, Mexico
§ National University of Colombia, Bogotá, Colombia
| University of Hildesheim, Hildesheim, Germany
Open Access


Insect population decline has been reported worldwide, including those of pollinators important for ecosystem services. Therefore, conservation actions which rely on available rigorous species distribution data are necessary to protect biodiversity. Niche modeling is an appropriate approach to distribution maps, but when it comes to bumble bees, few studies have been performed in South America. We modeled ecological niches of nine Colombian Bombus species with MAXENT 3.4 software using bioclimatic variables available from WorldClim. This resulted in maps for each species that show the potential distribution area at the present time. Modeled species maps accurately represent potential niches according to the description of bioclimatic conditions in the species’ habitat. We grouped the species into three clusters based on our results, as well as on distributional information from literature on the topic: High Mountain, Mid- Mountain and inter-Andean, and the Amazon and Eastern Plains Basin. Niche modeling depicted bumble bee species’ distribution in Colombia, the results of which can serve as a useful tool for conservation policies in the country.


biogeography, distribution, maxent, native bees, pollinators


Bumble bee species (Bombus spp.) are one of the widest studied taxa of bees in the world and Colombia regarding their ecological traits, the genetic composition of populations, pathogens, and their role in pollination services (Cure and Rodríguez 2007; Gamboa et al. 2015; Jaramillo-Silva et al. 2018; Lotta-Arevalo et al. 2020). Buzz pollination carried out by bumble bees allows for the pollination of several crops of economic importance in Colombia, such as Solanaceae like potatoes, tomatoes, paprika, etc. (Cure and Rodríguez 2007; Riaño et al. 2015; Vinícius-Silva et al. 2017) or Passifloraceae like passion fruit, gulupa, curuba, etc. (Pinilla-Gallego et al. 2016). Also, bumble bees are important for the pollination of many plant species in natural ecosystems like paramos and cloud forests, which are key for water regulation in Colombia (Rubio 2012). For example, they serve in the pollination of several species of “frailejon” of the genus Espeletia (Fagua and Bonilla Gómez 2005).

However, a drastic decline in insects has been documented worldwide (Hallmann et al. 2017; Seibold et al. 2019; Klink et al. 2020). A massive decline in the species distribution of most bumble bees has been reported, especially in developed regions like Europe and North America (Goulson et al. 2008). In Europe, there have been reports of declines since the 50’s. Just in the United Kingdom, three species out of twenty-five have gone extinct, and eight report major declines (Goulson 2003). Similar results have been reported in North America. In Illinois alone, four species have become locally extinct (Grixti et al. 2009). There is no population trend research of this kind performed in Colombia, and as one of the richest countries in terms of biodiversity (Rangel 2015), the country should increase research on its pollinators in order to develop preventive management of natural resources.

Biogeographic information about species is vital for conservation planning. The study of species distribution patterns has two complementary approaches: historical biogeography, which elucidates the causal mechanisms of current distribution, and ecological biogeography, which evaluates the ecogeographic factors that are currently shaping the distribution of species, looking for patterns in characteristics that are required for a species’ long-term survival. It looks at abiotic characteristics like humidity, temperature, or salinity and ecological characteristics like interactions with other organisms or genetical features. Research in this last approach is usually performed at a local scale a short time frame (Pérez-Malváez and Gutiérrez 2003).

Niche modeling is one of the methods used in ecological biogeography. It represents the fundamental niche without considering the realized niche, which requires detailed field research and confirmed species samples. Niche modeling aims to find the potential distribution area of a species (Soberon and Miller 2009). It is important to have good quality and quantity of information for the models, as they can be biased if the data for a species is scarce.

Thus, niche modeling can provide potential distribution maps for bumble bees in Colombia, updating and complementing the previous available maps, made by Abramovich et al. (2004). This information could be an important tool to use in both future studies with an evolutionary and ecological aim, as well as in conservation efforts by informing environmental conservation plans, political decisions, and public policies in Colombia, and by strengthening existing efforts like the Colombian Pollinators Initiative (CPI) (Nates-Parra 2016). As the importance of bee conservation continues to escalate, conservation actions are necessary to protect bees’ biodiversity (Sárospataki et al. 2005).

Materials and methods

Study area

The geographical extent used for each species in the modeling process covers the entire continental Colombian territory, from the north (13°23.73'N) to the south (4°13.75'N) and from the west (81°44.13'W) to the east (66°50.63'W). No buffer was used.

There are nine species from the genus Bombus reported for Colombia: B. pauloensis Friese, 1913 (formerly B. atratus Franklin, 1913), B. excellens Smith, 1879, B. funebris Smith, 1854, B. hortulanus Friese, 1904, B. melaleucus Handlirschi, 1888, B. pullatus Franklin, 1913, B. robustus Smith, 1854, B. rubicundus Smith, 1854, and B. transversalis Oliver, 1789. Across the paramo ecosystem, four species can be found with different altitudinal distributions: B. funebris between 2500 and 4750 m, B. hortulanus between 2100 and 3600 m, B. robustus between 2100 and 3800 m, and B. rubicundus between 2500 and 3900 m (Pinilla-Gallego et al. 2016). Along the sub-Andean forest or the cloud forest, two species can be found: B. excellens between 1500 and 2600 m and B. melaleucus between 450 and 2100 m. The most abundant species in the low mountain strata is B. pauloensis. However, it has a wide altitudinal range between 150 and 3600 m (Liévano et al. 1991). Finally, two different species can be found across tropical forests with warm and wet climates: B. pullatus between 120 and 3500 m (especially in the foothills) and B. transversalis between 180 and 1100 m (restricted to the amazon forest and its foothills).

Occurrence data

We obtained occurrence points from the Wild Bee Research Lab of the Universidad Nacional de Colombia, Bogotá (LABUN), the Insect Collection of Universidad del Quindío (CIUQ), the Francisco Luis Gallego Entomological Museum (MEFLG), Universidad Nacional de Colombia, Medellín, the Bee Collection of the Universidad Militar Nueva Granada, the Entomological Museum of the School of Agronomy at Universidad Nacional de Colombia (UNAB), and the database from the Global Biodiversity Information Facility (GBIF). Occurrence points were adjusted using the altitude reference reported by Pinilla-Gallego et al. (2016). All occurrence point data were originally collected from 1938 until 2020, throughout different seasons and derived from various methods (See Suppl. material 1). The senior author checked all occurrence records to prevent recognizable errors in georeferencing and taxonomy. The number of points used for each specie was: B. excellens 46; B. funebris 150; B. hortulanus 1130; B. melaleucus 47; B. pauloensis 2128; B. pullatus 549; B. robustus 243; B. rubicundus 426; B. transversalis 43 (Fig. 1).

Figure 1. 

Distribution of the points used in the modeling of each specie.

Environmental variables

We used 19 environmental data layers for modeling (Table 1), which were downloaded from the WorldClim V. 2.1 website (Fick and Hijmans 2017). The layers represent the climate average from the year 1970 to the year 2000 at a spatial resolution of 30 Arc seconds (0.93 km2).

Table 1.

Bioclimatic variables included in modeling were obtained from WorldClim.

Temperature Precipitation
1 Annual Mean Temperature 12 Annual Precipitation
2 Mean Diurnal Range (Mean of monthly (max. temp.-min. temp.)) 13 Precipitation of Wettest Month
3 Isothermality (BIO2/BIO7) (×100) 14 Precipitation of Driest Month
4 Temperature Seasonality (standard deviation ×100) 15 Precipitation Seasonality (Coefficient of Variation)
5 Max. temperature of Warmest Month 16 Precipitation of Wettest Quarter
6 Min. Temperature of Coldest Month 17 Precipitation of Driest Quarter
7 Temperature Annual Range (5–6) 18 Precipitation of Warmest Quarter
8 Mean Temperature of Wettest Quarter 19 Precipitation of Coldest Quarter
9 Mean Temperature of Driest Quarter
10 Mean Temperature of Warmest Quarter
11 Mean Temperature of Coldest Quarter

Niche models

QGIS 2.8 (QGIS Development Team 2016) was used to clip and convert environmental layers, enabling their use in MAXENT. Niche models were developed using MAXENT 3.4 (Phillips et al. 2020). Specific settings were set at the following defaults: a maximum of 500 iterations, 10% test point (CrossValidate), without extrapolation and cumulative, keeping a limit convergence of 0.00001 and prevalence of 0.5, maximum number of background points 10000.

Two different statistical analyses were performed per species: the Jackknife test to evaluate the weight or importance of each variable (Timaná de la Flor and Romero 2015; Phillips and AT&T Research 2017), and the AUC (Area Under the Curve) test. Jackknife test is a method for validating the samples and model. It shows the representative variables for modeling each specie (Shcheglovitova and Anderson 2013). According to the results of the Jackknife test, we selected the three variables with the highest percentage of importance, whose values depend on each species. To obtain the AUC, the specific settings given above were chosen according to established knowledge about bumble bees’ altitudinal distribution (Pinilla-Gallego et al. 2016). Subsequently, an ROC (Receiver Operating Characteristics) curve was created, and the AUC was calculated. The AUC value reflects the model’s accuracy or its capacity for prediction. The AUC value moves between 0 and 1, with 1 representing a perfect prediction. Values over 0,9 are considered strong (Pliscoff and Fuentes-Castillo 2011). Finally, ArcMap 10.5 (ESRI 2015) was used for reclassifying and visually representing the ecological niche modeling results.


The models showed good results for all nine bumble bee species, with AUC values above 0.9, although the number of occurrences for Bombus excellens, Bombus melaleucus, and Bombus transversalis was low. Bombus pauloensis and Bombus hortulanus had the highest number of occurrences, being common in collections (Table 2). Based on our results and distributional information from scientific literature, we divided the bumble bee species into three distributional groups: high mountain bumble bees, mid-mountain/inter-andean bumble bees, and Amazon and eastern Plains basin bumble bees.

Table 2.

Bumble bee species data for the area under the curve (AUC obtained by niche modeling), number of occurrences used for modeling, and the distributional group selected for each species according to our results and previous results by Pinilla-Gallego et al. (2016).

Species AUC 1970-2000 Distributional group
Bombus excellens 0.951 Mid-mountain and inter-Andean
Bombus funebris 0.976 High Mountain
Bombus hortulanus 0.977 High Mountain
Bombus melaleucus 0.965 Mid-mountain and Inter-Andean
Bombus pauloensis 0.911 Mid-mountain and Inter-Andean
Bombus pullatus 0.917 Mid-mountain and Inter-Andean
Bombus robustus 0.983 High Mountain
Bombus rubicundus 0.973 High Mountain
Bombus transversalis 0.999 Amazon and Eastern Plains Basin

High Mountain bumble bees

The potential distribution areas for B. funebris, B. hortulanus, B. robustus, and B. rubicundus occurred only in high mountain departments, along a range of different small and fragmented areas, with a high probability of occurrence in the central part of the Eastern Andes Range (Fig. 2). B. funebris represented the most restricted and fragmented species. B. hortulanus occurrences had the lowest altitudinal points for this group. B. robustus data was similar to B. hortulanus but included some areas in the northern part of the Eastern and the Central Ranges. B. rubicundus showed similar potential distribution to B. robustus but with more restricted and fragmented areas.

Figure 2. 

Potential distribution maps for High Mountain bumble bees A B. funebris B B. hortulanus C B. robustus D B. rubicundus.

Mid-Mountain and Inter-Andean bumble bees

This group of species showed a wider altitudinal and potential distribution than the high mountain species, along and between the three mountain ranges (Fig. 3). B. excellens showed a higher probability along most of the Central and Eastern mountain ranges, as well as in the mid-mountain area of the Sierra Nevada de Santa Marta.

Figure 3. 

Potential distribution maps for Mid-Mountain and Inter-Andean bumble bees. A B. excellens B B melaleucus C B. pauloensis D B. pullatus.

For B. melaleucus, its potential distribution was indicated along the three mountain ranges. B. pauloensis exhibited a continuous distribution over the Andean region. B. pullatus was concentrated along low altitudinal areas of the inter-Andean valleys.

Amazon and Eastern Plains Basin bumble bees

Bombus transversalis is the only species with a potential distribution along the Amazon Forest, the Andean foothills, and the Orinoquia region. The eastern part of Meta presented the highest altitudinal points for the species at 1150 m. The northern area of the biogeographic region of Chocó showed optimal bioclimatic conditions for B. transversalis, but no occurrence point was found or used there for this model (Fig. 4).

Figure 4. 

Potential distribution map for Amazon and Eastern Plains Basin bumble bees, B. transversalis.

The variables with the highest weight for each species are presented in Table 3.

Table 3.

Bioclimatic variables had the highest weight in the model for each species, according to the results obtained by the Jackknife test.

Species (Distributional group) Variables with the highest weight (Percent of contribution)
1st 2nd 3rd
Bombus excellens (Mid-Mountain and inter-Andean) 8. Mean Temperature of Wettest Quarter 94.35% 1. Annual Mean Temperature 94.27% 6. Min. Temperature of Coldest Month 94.04%
Bombus funebris (High Mountain) 4. Temperature Seasonality 95.50% 6. Min. Temperature of Coldest Month 94.54% 1. Annual Mean Temperature 94.03%
Bombus hortulanus (High Mountain) 11. Mean Temperature of Coldest Quarter 96.72% 1. Annual Mean Temperature 96.63% 10. Mean Temperature of Warmest Quarter 96.61%
Bombus melaleucus (Mid-Mountain and inter-Andean) 5. Max. Temperature of Warmest Month 95.20% 11. Mean Temperature of Coldest Quarter 95.16% 1. Annual Mean Temperature 95.14%
Bombus pauloensis (High Mountain) 6. Min. Temperature of Coldest Month 91.30% 11. Mean Temperature of Coldest Quarter 91.18% 8. Mean Temperature of Wettest Quarter 91.09%
Bombus pullatus (Mid-Mountain and inter-Andean) 11. Mean Temperature of Coldest Quarter 83.91% 8. Mean Temperature of Wettest Quarter 83.53% 4. Temperature Seasonality 76.67%
Bombus robustus (High Mountain) 5. Max. Temperature of Warmest Month 93.74% 6. Min. Temperature of Coldest Month 93.72% 8. Mean Temperature of Wettest Quarter 93.51%
Bombus rubicundus (High Mountain) 6. Min. Temperature of Coldest Month 96.70% 1. Annual Mean Temperature 96.70% 8. Mean Temperature of Wettest Quarter 96.48%
Bombus transversalis (Amazon and Eastern Plains Basin) 13. Precipitation of Wettest Month 98.88% 11. Mean Temperature of Coldest Quarter 97.87% 6. Min. Temperature of Coldest Month 96.94%


These potential distribution maps for bumble bees in Colombia improve the previous maps available (Abrahamovich et al. 2004) for the species, with more precise areas of potential distribution and a visual schema for patterns already described in the literature (Pinilla-Gallego et al. 2016).

The four high-mountain species are associated with the paramo and high Andean ecosystems (Pinilla-Gallego et al. 2016). Paramo ecosystems are characterized as biogeographical islands, highly isolated areas with a great deal of endemism and low genetic flow between populations (Lotta-Arevalo et al. 2020). As a result, the potential distribution of these species is highly fragmented. B. excellens, B. melaleucus, and B. pullatus are associated with mountain forests, and it is remarkable that, even though they showed a wide suitable area (Fig. 3A, B, D, respectively), there is a low number of occurrences. This indicates a small population size and high vulnerability. For example, B. excellens can only be found in the cloud forest, an ecosystem with an accelerated rate of deforestation (Nates-Parra 2006). B. pauloensis (Fig. 3C) shows a wide altitudinal and potential distribution consistent with its high plasticity in the habitat chosen for nesting (Liévano et al. 1991; Nates-Parra et al. 2006), and this is reflected by its wide distribution along the mountain areas. Its potential distribution will help conservation planning, as B. pauloensis positively impacts national agricultural productivity by pollinating staple foods such as fruits and vegetables (Cure and Rodríguez 2007; Riaño et al. 2015; Poveda et al. 2018) due to its adaptability to disturbed habitats. According to its biogeographical history, for B. transversalis the Andes Mountain Range represents a physical barrier for emigration. Thus, the potential distribution area in the biogeographical region of Chocó is not likely (Abrahamovich et al. 2004). The eastern area of Meta stands out as an expansion of previously recorded areas to include the low areas of the Andean foothills.

The most important bioclimatic variables for bumble bees in Colombia are related to temperature (Table 3). Thus, drastic temperature changes put bumble bees in a condition vulnerable to threat under a climate change scenario (Gonzalez et al. 2021). Also, the potential distribution for most species overlaps with densely populated areas. Therefore, these maps are a tool for the detailed location of bumblebee biodiversity for conservation planning as bumble bee’ populations and their associated services will probably soon be reduced (Williams and Osborne 2009; Pinilla-Gallego et al. 2016).


The potential distribution for Colombian bumble bee species is reported here. Seven of them show a restricted distribution, shaped mainly by temperature restrictions. These results are a direct contribution to knowledge about Colombian bumble bee species, constituting useful knowledge for conservation, territorial planning, protection plans, and environmental management. Thus, by obtaining the most suitable areas for a species, this study can provide the ideal location to breed and reproduce the native species. These species contribute to improving the agricultural productivity of several regions by pollination. Likewise, this information can help avoid the introduction of foreign species for this purpose (Nates-Parra 2016). It is necessary to increase research efforts on bees for them to be included in conservation planning.


We would like to thank M. Sc. Andrea Lorena García and M. Sc. Daniela Hoyos-Benjumea from the Insect collection at Universidad del Quindío, Ph.D. Diego Riaño from the bee collection in Universidad Militar Nueva Granada, Ph.D. Francisco Serna from the entomological museum Universidad Nacional, School of Agronomy, Bogotá, and Ph.D. Sergio Orduz and biol. Carlos Londoño from the entomological museum Francisco Luis Gallego for their contributions to occurrence data. We are indebted to Ph.D. Sydney A. Cameron for comments and suggestions that improved this manuscript.

We acknowledge financial support by Stiftung Universität Hildesheim.


  • Abrahamovich AH, Díaz NB, Morrone JJ (2004) Distributional patterns of the Neotropical and Andean species of the genus Bombus (Hymenoptera: Apidae). Acta zoológica mexicana 20: 99–117.
  • Cure JR, Rodríguez D (2007) Efecto de Bombus atratus (Hymenoptera: Apidae) sobre la productividad de tomate (Lycopersicon esculentum Mill.) bajo invernadero en la Sabana de Bogotá, Colombia. Agronomía Colombiana 25: 62–72.
  • ESRI (2015) ArcGIS ArcMap 10.5.
  • Fagua JC, Bonilla Gómez MA (2005) Ecología de la polinización de Espeletia grandiflora en el Parque Chingaza. In: Estrategias adaptativas de plantas del páramo y del bosque altoandino en la Cordillera Oriental de Colombia. Universidad Nacional de Colombia.
  • Fick SE, Hijmans RJ (2017) WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37: 4302–4315.
  • Gamboa V, Ravoet J, Brunain M, Smagghe G, Meeus I, Figueroa J, Riaño D, de Graaf DC (2015) Bee pathogens found in Bombus atratus from Colombia: A case study. Journal of Invertebrate Pathology 129: 36–39.
  • Gonzalez V, Cobos ME, Jaramillo-Silva J, Ospina-Torres R (2021) Climate change will reduce the potential distribution ranges of Colombia’s most valuable pollinators. Perspectives in Ecology and Conservation 19(2): 195–206.
  • Goulson D (2003) Bumblebees: Their Behaviour and Ecology. Oxford University Press, 260 pp.
  • Jaramillo-Silva J, Ospina-Torres R, Henao-Cardenas MM, Rodriguez A, Parra G (2018) Passion fruit in Colombia. In: The pollination of cultivated plants. A compendium for practitioners. Balboa, Ancon, Republic of Panama, 226–232.
  • Liévano A, Ospina-Torres R, Nates-Parra G (1991) Distribución altitudinal del género Bombus en Colombia (Hymenoptera: Apidae). Trianea 4: 541–550.
  • Lotta-Arevalo I, Vargas-Ramírez M, Parra G, Matta N, Ospina-Torres R (2020) Accediendo al pasado: uso de especímenes de colección como fuentes de información genética para el género Bombus (Hymenoptera: Apidae). Revista de Biología tropical 68: 394–414.
  • Nates-Parra G [Ed.] (2016) Iniciativa colombiana de polinizadores: Abejas ICPA. In: Iniciativa colombiana de polinizadores. Departamento de Biología, Universidad Nacional de Colombia, Bogotá D. C., Colombia, 42 pp. http://localhost:8080/handle/11438/8800 [April 16, 2021]
  • Nates-Parra G, Parra A, Rodríguez A, Baquero P, Vélez D (2006) Abejas silvestres (Hymenoptera: Apoidea) en ecosistemas urbanos: Estudio en la ciudad de Bogotá y sus alrededores. 1: 77–84.
  • Pérez-Malváez C, Gutiérrez RR (2003) Las ideas biogeográficas y su presencia en una revista mexicana: La Naturaleza. Llull: Revista de la Sociedad Española de Historia de las Ciencias y de las Técnicas 26: 207–244.
  • Phillips SJ (2017) A brief tutorial on Maxent., 1–39.
  • Pliscoff P, Fuentes-Castillo T (2011) Modelación de la distribución de especies y ecosistemas en el tiempo y en el espacio: una revisión de las nuevas herramientas y enfoques disponibles. Revista de Geografía Norte Grande 48: 61–79.
  • Poveda CA, Riaño Jimenez D, Aguilar Benavides L, Cure JR (2018) Eficiencia de polinización de colonias huerfanas de Bombus atratus (Hymenoptera: Apidae) en fresa (Fragaria × ananassa) bajo cubierta. Acta Biológica Colombiana 23: 73–79.
  • QGIS Development Team (2016) QGIS Geographic information system 2.8.9. Open Source Geospatial Foundation Project.
  • Rangel O (2015) La biodiversidad de Colombia: significado y distribución regional. Revista de la Academia Colombiana de Ciencias Exactas, Físicas y Naturales 39: 176–200.
  • Riaño D, Pacateque J, Cure JR, Rodríguez D (2015) Pollination behavior and efficiency of Bombus atratus Franklin in sweet peppers (Capsicum annum L.) grown in a greenhouse. Revista Colombiana de Ciencias Hortícolas 9: 259–267.
  • Rubio Fernández D (2012) Disponibilidad, uso y preferencia por los recursos florales en una comunidad de abejorros (Hymenoptera: Apidae: Bombus) en el páramo de Chingaza. Universidad Nacional de Colombia. [March 16, 2021]
  • Vinícius-Silva R, Parma D de F, Tostes RB, Arruda VM, Werneck M do V (2017) Importance of bees in pollination of Solanum lycopersicum L. (Solanaceae) in open-field of the Southeast of Minas Gerais State, Brazil. Instituto de Botânica. [November 12, 2022]

Supplementary material

Supplementary material 1 

Points used for the Niche modeling of Bumblebee species (Hymenoptera, Apidae, Bombus) in Colombia

Laura Rojas-Arias, Daniel Gómez-Morales, Stephanie Stiegel, Rodulfo Ospina-Torres

Data type: COL (excel document).

Explanation note: Ocurrence data used for the modeling.

This dataset is made available under the Open Database License ( The Open Database License (ODbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.
Download file (120.09 kb)
login to comment