Climate‐induced distribution dynamics of Plebeia flavocincta, a stingless bee from Brazilian tropical dry forests

Abstract Aim The objective of this study is to estimate the current potential geographic distribution of Plebeia flavocincta and to evaluate the influence of climate on the dynamics of suitable habitat availability in the past and in the future. Location Northeast region of Brazil and dry forest areas. Methods The habitat suitability modeling was based on two algorithms, two global circulation models, and six different scenarios. We used this tool to estimate the areas of occurrence in the past (Last Interglacial and Last Glacial Maximum), in the present, and in the future (years 2050 and 2070). Results According to the models, P. flavocincta had great dynamics in the availability of suitable habitats with periods of retraction and expansion of these areas in the past. Our results suggest that this taxon may benefit in terms of climate suitability gain in Northeast Brazil in the future. In addition, we identified high‐altitude areas and the eastern coast as climatically stable. Conclusion The information provided can be used by decision makers to support actions toward protecting and sustainably managing this taxon. Protection measures for this taxon are particularly important because this insect contributes to the local flora and, although our results indicate that the climate may favor this taxon, other factors can negatively affect it, such as high levels of habitat loss due to anthropogenic activities.

of Northeast Brazil (Figure 1). This taxon occurs in tropical dry forests (Zanella & Martins, 2003) and interacts with approximately 60 plant species (Costa, 2014), and can be considered an important visitor of the native flora.
Tropical dry forests of northeastern Brazil, also known as Caatinga, consist of an ecoregion composed predominantly of xerophilic, woody, thorny, and deciduous vegetation structures (Costa, de Araújo, & Lima-Verde, 2007;Sampaio, 1995), and a hot and dry climate (Veloso, Rangel-Filho, & Lima, 1991), extending over more than 750,000 km 2 (Silva, Barbosa, Leal, & Tabarelli, 2017). The unique climatic conditions of the Caatinga are responsible for the existence of unique biodiversity adaptations, which allow for a large number of species adapted to heat and drought. Caatinga areas are undergoing extensive processes of environmental alteration and deterioration caused by unsustainable resource use (Leal, Tabarelli, & da Silva, 2003;Silva et al., 2017). In addition, there is a high risk of desertification in the Caatinga, resulting from several factors, including human activities and climatic variations (Santana, 2007). Different initiatives have projected future climate change scenarios, with most predicting an increase in temperature and a decrease in precipitation in this region (Marengo, Alves, Beserra, & Lacerda, 2011;Torres, Lapola, & Gamarra, 2017). According to the 5th assessment report of the Intergovernmental Panel on Climate Change (AR5 IPCC 2014), climate scenarios suggest that the main threat to the Caatinga is increasing aridity, which could lead to an intensification and prolongation of the dry seasons. However, at least for endemic birds and reptiles, most species analyzed using habitat suitability models (HSMs) were projected to gain climatically suitable areas in remnants of natural Caatinga vegetation (Oliveira, Araújo, Rangel, Alagador, & Diniz-Filho, 2012).
There is extensive evidence of a decline in pollinators around the world, with habitat fragmentation/loss and climate change as the main causes of this decline (Potts et al., 2010). Climate change can affect bees on different levels, such as behavior, physiology, quality of the floral environment, population dynamics, and the emergence of new interactions between species and parasites and pathogens (Le Conte & Navajas, 2008). Small bee species are more likely to be affected by changes in climatic conditions because they lose heat faster than larger bees due to their small body size (Oyen & Dillon, 2018). This is supported by Costa (2014) and Silva, Meneses, and Freitas (2019), who observed that P. flavocincta prefers to start its flight activity during periods of higher temperatures.
The spatial distribution of P. flavocincta is not precisely known, which hinders the development of conservation and management plans for this taxon. Habitat suitability models have been used to both evaluate the niche dynamics of a species and determine the most important environmental variables that affect its potential distribution (Pearman, Guisan, Broennimann, & Randin, 2008). The applications of this tool for conservation purposes are based on the projection of climatically suitable areas in the present and in the future or on the investigation of the potential historical distribution of a species. This tool is increasingly used to predict the effects of climate change over time on species distribution (Franklin, 2010).
The objective of this study is to estimate the current potential geographic distribution of P. flavocincta and evaluate the influence of climate on the dynamics of suitable habitat availability in the past and in the future. We focused on occurrence records of what is now F I G U R E 1 Workers of Plebeia flavocincta at the entrance of the nest in a rational wooden hive F I G U R E 2 Study area used for model projections. The thick black line delimits the Northeast region of Brazil; the green area delimits the Caatinga domain; the gray gradient represents altitude-darker colors, higher altitude; rivers are in blue; and green triangles are the occurrence records of Plebeia flavocincta used for this study taxonomically recognized as P. flavocincta and asked how climatic fluctuations could affect the distribution dynamics of the taxon in the past and in the future.
The Caatinga constitutes most of this area (Figure 2), which has a semiarid climate with precipitation of 500-800 mm/year (Velloso, Sampaio, & Pareyn, 2002). Precipitation is strongly seasonal (during only 3-4 months of the year, between December and March) and has a high interannual variability, which leads to severe recurrent and irregular droughts (Brito et al., 2018).

| Georeferenced occurrences
We gathered data from the digital platforms Global Biodiversity Information Facility (GBIF; www.gbif.org, accessed in February 2019), SPECIESLINK (www.splink.org.br, accessed in February 2019), and complementary data collected in the field, totaling 615 records.
We verified the information about who identified each specimen in each record found in the digital databases, and we kept only entries made by specialists. We also applied a filter to remove duplicates and nongeoreferenced data. Thus, our final database contained 52 occurrence records (30 from digital database and 22 from our field collections) for P. flavocincta in seven Brazilian states ( Figure 2). Most of the points are located in Caatinga areas.

| Climate data
The species distribution models were calibrated using current climate data from the WorldClim database (www.world clim.org, accessed in February 2019) with a resolution of ~1 km 2 . We tested which of the 20 variables were most strongly correlated with each other (r > .75) through the "removeCollinearity" function of the virtualspecies package in R (Leroy, Meynard, Bellard, & Courchamp, 2016;R Development Core Team, 2019) and selected nine bioclimatic variables: mean diurnal range (bio 2), isothermality (bio 3), temperature seasonality (bio 4), maximum temperature of the warmest month (bio 5), annual precipitation (bio 12), precipitation seasonality (bio 15), precipitation of the warmest quarter (bio 18), precipitation of the coldest quarter (bio 19), and altitude.

| Habitat suitability model
The HSM was based on records of presence and pseudoabsence. We generated three sets of pseudoabsence data containing ten times the number of randomly distributed presence data points (Chefaoui & Lobo, 2008). We used two algorithms: generalized linear model (GLM) (McCullagh & Nelder, 1989) and maximum entropy (MAXENT) (Phillips, Anderson, & Schapire, 2006), which have been commonly used for generating robust responses (Li & Wang, 2013). The databases were partitioned, with 80% of the original data used for model calibration and 20% for evaluation. This procedure was repeated 10 times to ensure that the predictive accuracy of the model was not affected by the random partition strategy. In addition, we evaluated the quality of these runs according to two different metrics: true skill statistics (TSS, Allouche, Tsoar, & Kadmon, 2006) and receiver operating characteristic curve (ROC; Phillips et al., 2006).

| Model performance and environmental variables
In general, the different runs of the models exhibited high evaluation scores, that is, values that reflect the quality of the models (TSS and ROC > 0.8). The mean diurnal range (bio 2) and annual precipitation (bio 12) were the most important environmental variables for HSMs ( Figure 3). According to the models, in the LIG (~120 kya), the eastern coast of the Northeast region was relatively colder and drier compared to the current mean diurnal temperature range and annual precipitation. For the LGM (~21 kya), this same area had higher mean temperature ranges compared to the present; in addition, the central Caatinga had higher annual precipitation rates. annual precipitation (bio 12), Northeast Brazil will become even drier compared to today, as it was at 120 kya, mainly along the entire eastern coast, where the precipitation rate may decrease >600 mm.

| Species distribution over time
According to our models, the distribution of suitable habitats for P. flavocincta during the LIG was greater than the current potential

| Stability of the areas
Consensus maps of all scenarios were overlapped to determine common areas that represented potential occurrence regardless of the scenario (which we considered climatically stable and, therefore, we call refuge areas; Figure 6). When analyzed with all models (past, present, and future, totaling 24 models), the refuge areas are found in a low region of the São Francisco River Valley, between the states of Sergipe and Alagoas, on the eastern coast (in the states of Pernambuco, Alagoas, and Sergipe) and a small transition area of Atlantic Forest and Caatinga of Paraíba (Figure 6a). When only the projection models for the future were analyzed, several climatically suitable areas were also observed in low areas of the São Francisco River Valley, in areas of the Serra da Ibiapaba, Serra da Mata, and Maciço de Baturité mountain ranges (state of Ceará), transition area of Atlantic Forest and Caatinga of Paraíba and along the entire eastern coast, in several states of the Northeast region (Figure 6b).

| D ISCUSS I ON
The distribution dynamics of P. flavocincta in the past were intense, as the total area of climatically suitable habitats projected for the oldest period (LIG, ~120,000 years ago) was much larger than the area of climatically suitable habitats projected for the present. In addition to being the environmental variables that most contributed to explaining the models, the mean diurnal range (bio 2) and annual precipitation (bio 12) varied greatly over time. The projected locations of suitable habitat for P. flavocincta seem to be associated with higher temperatures and lower precipitation rates. Some studies show that bee flight activity may be limited depending on temperature and precipitation rate (Hilário, Ribeiro, & Imperatriz-Fonseca, 2012;Silva et al., 2019). Due to a small body size, workers lose heat faster and therefore need higher temperatures for foraging activities (Oyen & Dillon, 2018). Researchers  Flight activity is directly related to the dispersal ability of a species. Our models projected for the past suggest an intense dynamic in the distribution of suitable habitats with signs of retraction-expansion through which populations potentially became isolated and that may have left genetic signatures. According to paleodistribution and palynological validation models, the potential historical distribution of Caatinga areas was more restricted during the LGM (Werneck, Costa, Colli, Prado, & Sites, 2011). This retraction phase of the Caatinga during the supposedly dry and cold period of the LGM was followed by an expansion with small fluctuations until reaching the current area. This phenomenon may explain the distribution dynamics of P. flavocincta during periods of vegetation retraction and expansion. In fact, it is possible to observe variations in P. flavocincta populations regarding behavior, worker coloration, and hive morphology (personal observations) that may be associated with this dynamic. However, for another stingless bee species endemic to the region, Partamona rustica Pedro & Camargo, 2003, its potential distribution has been continuously expanding since the LIG (Miranda et al., 2016). For this reason, more specific studies on population variability (i.e., wing morphometry and genetic analyses) are needed to confirm these assumptions.
Our results also suggest that P. flavocincta will find a larger area of climatically suitable habitat in the 2050s and especially in the Due to the lower number of samples in the central Caatinga region, these results may underestimate the potential distribution of P. flavocincta in this area. In turn, this region is known to be susceptible to desertification (Salazar, Nobre, & Oyama, 2007) due to the joint effects of increasing temperatures, reduced precipitation and soil degradation (Darhoh, 1998;Geist & Lambin, 2004;Sivakumar, 2007).
From a climate change perspective, the area size of a domain may change (Salazar et al., 2007), and new climatically suitable areas may emerge, causing P. flavocincta to further expand its area of occurrence. However, underlying the spatial structure of climatic impacts is the pressures of human occupation, land degradation, and biological interactions that may prevent the taxon from occupying these areas. Despite covering an area of approximately 10% of the national territory, only 14% of the Caatinga is protected by law as national conservation areas (Gariglio, Sampaio, Cestaro, & Kageyama, 2010;Sparovek, Berndes, Klug, & Barretto, 2010); that is, its biodiversity has been largely ignored by conservation policies (Leal, Silva, Tab arelli, & Lacher, 2005). Many of the new climatically suitable areas for P. flavocincta (2050 and 2070) will coincide with Conservation Units (www.mma.gov.br/areas -prote gidas/ cadas tro-nacio nal-de-ucs/map). These conservation areas are important to contain or reduce deforestation and allow the taxon to shelter in protected areas.
F I G U R E 6 Areas of potential climatic stability resulting from the overlapping of suitability maps for Plebeia flavocincta. The areas in yellow represent locations indicated by all scenarios/models as suitable for the species. (a) past, present, and future, "All models" = 24; (b) present and future only, "All models" = 18 The refuge areas estimated by our results are located in the mountains of Ceará, in the central portion of Paraíba, on the eastern coast, and in the lower region of the São Francisco River Valley, among the states of Sergipe, Alagoas, and Pernambuco. These refuge areas are sites where climatic changes are less intense, remaining suitable habitats for the taxon over time and, therefore, possibly representing priority conservation areas. In addition to delimiting priority areas, the planting of native species that provide food resources for P. flavocincta, such as Mimosa tenuiflora (Willd.) Poir. (Fabaceae), may represent an important step for species conservation (Costa, 2014). Such plantings may be combined with the recovery of degraded areas and the creation of ecological corridors, especially in times of climatic change and scarcity of floral resources in the dry season. Recently, das Chagas, Lucas, and de Almeida Vieira (2020) projected a potential reduction in the climatically suitable areas for M. tenuiflora when comparing the current and the future (2070) scenarios. The reduction of climatically suitable areas was discussed as potentially causing spatial incompatibility in pollination interactions (Gérard, Vanderplanck, Wood, & Michez, 2020).
However, little is known about climate change impact on plant-bee interactions (but see Giannini et al., 2013); moreover, species' ability to adapt and persist is often neglected. It is also necessary to consider that most species of the tribe Meliponini nest in the hollows of living trees (Nogueira-Neto, 1997) and are dependent on forest to nidify (Giannini et al., 2015(Giannini et al., , 2017. In areas of Caatinga, P. flavocincta prefers to nest in Cenostigma pyramidale (Tul.) Gagnon & GP Lewis (Fabaceae), a native tree species popularly known as catingueira and widely used for nesting by other stingless bee species (Martins, Laurino, Koedam, & Fonseca, 2004).

| CON CLUS ION
Our results contribute to the interpretation of biogeographic scenarios in Caatinga areas and reinforce the need for management and conservation plans for P. flavocincta in priority areas, where the taxon potentially finds stable climatic conditions in past and future scenarios. Protection measures for this taxon are particularly important because this insect contributes to the local flora as it exhibits generalist behavior. Moreover, P. flavocincta contributes to local communities as a source of income from honey production, as it is the second most farmed bee species in the state of Rio Grande do Norte (Maia, Jaffe, Carvalho, & Fonseca, 2015) and the third most farmed in the state of Ceará (Felix, 2015).

ACK N OWLED G M ENTS
We thank the Fundação de Amparo à Pesquisa do Estado do Pará-FAPESPA (PhD scholarship), and the Coordenação de

CO N FLI C T O F I NTE R E S T S
The authors have declared that no competing interests exist.

DATA AVA I L A B I L I T Y S TAT E M E N T
All relevant data are within the manuscript and its Supporting Information files (Appendix S1).