Different Responses in Geographic Range Shifts and Increase of Niche Overlap in Future Climate Scenario of the Subspecies of Melipona quadrifasciata Lepeletier

Pollination is considered a key element of global biodiversity conservation, because it provides a crucial ecosystem service, as without pollination, many plants would not be able to reproduce (Klein et al., 2007). Ollerton et al. (2011) estimated that approximately 87.5% of world’s flowering plants require animal pollination for fruit and seed production. Pollination is also considered fundamental for food production because one-third of world’s leading food crops depends on animal pollination (Klein et al., 2007). In Brazil, Giannini et al. (20015) demonstrated that about 89% of the crop species are pollinated by bees. Abstract Climate change is suggested to be one of the possible drivers of decline in pollinators. In this paper, we applied an ecological niche model to modeling distributional responses in face of climate changes for the subspecies of Melipona quadrifasciata Lepeletier. This species is divided into two subspecies based on difference in the yellow tergal stripes, which are continuous in M. q. quadrifasciata and interrupted in M. q. anthidioides. The geographic distribution of each subspecies is also distinct. M. q. quadrifasciata is found in colder regions in the Southern states of Brazil, whereas M. q. anthidioides is found in habitats with higher temperatures, suggesting that ecological features, such as adaption to distinct climatic conditions may take place. Thus, the possibility of having different responses in geographic range shifts to future climate scenario would be expected. This study aimed to investigate the effects of climate changes on the distribution of the two M. quadrifasciata subspecies in Brazil, using an ecological niche model by the MaxEnt algorithm. Our results indicate that the subspecies showed clear differences in geographic shift patterns and increased climate niche overlap in the future scenarios. M. q. anthidioides will have the potential for an increase of suitable climatic conditions in the Atlantic Forest, and towards the Pampa biome, while M. q. quadrifasciata will suffer a reduction of adequate habitats in almost all of its current geographic distribution. Given the potential adverse effects of climate changes for this subspecies, conservation actions are urgently needed to avoid that it goes extinct. Sociobiology An international journal on social insects


Introduction
Pollination is considered a key element of global biodiversity conservation, because it provides a crucial ecosystem service, as without pollination, many plants would not be able to reproduce (Klein et al., 2007).Ollerton et al. (2011) estimated that approximately 87.5% of world's flowering plants require animal pollination for fruit and seed production.Pollination is also considered fundamental for food production because one-third of world's leading food crops depends on animal pollination (Klein et al., 2007).In Brazil, Giannini et al. (20015) demonstrated that about 89% of the crop species are pollinated by bees.
There is evidence of a global decline in pollinators, mainly honey bees (Becher et al., 2013), which points to a more important role of wild bees as providers of this ecological service.Wild bees, especially social stingless bees, are pollinators of many native plant species (Imperatriz-Fonseca et al., 2012) and effective crop pollinators (Klein et al., 2007).
Among the stingless bees, the Melipona genus has the largest number of species (Moure et al., 2007), and can be found throughout the Neotropical region, from Mexico to Misiones, Argentina, with highest diversity in the Amazon basin (Silveira et al., 2002).So far, declines of Melipona species in Brazil have gone unnoticed, until three species (Melipona bicolor schenkii Gribodo, Melipona marginata obscurior Moure, Melipona quadrifasciata quadrifasciata Lepeletier) were reported to be threatened in the Brazilian state of Rio Grande do Sul (Fontana et al., 2003).Currently, there are eight Melipona species reported to be threatened at national and/or state level (Tossulino et al., 2006;ICMBio, 2016).
Multiple anthropogenic pressures, including shortage of food sources and nesting sites due to habitat loss, aggressive agricultural practices and the spread of the alien species Apis mellifera L., are the main factors responsible for the observed population declines (Vanbergen et al., 2013).Climate change also has been considered one of the possible drivers of the pollinators decline, because climatic variables provide general conditions for their occurrence and performance, according to their physiological limits (Dupont et al., 2011).Mainly temperature and precipitation have been shown to play important roles in plant-pollinator relationships and may provide critical insights to help explain the potential effects of climate change on these interactions (Devoto et al., 2009).Shifts in phenology induced by climate change have the potential to disrupt the temporal overlap between pollinators and their floral food resources.As pollinators will experiencing loss of some of their food plants, they are likely to suffer population declines, as relying on fewer plant species will expose them to lower overall densities of flowers and greater temporal and spatial variation in food supply (Memmott et al., 2007).
There are three general expectations for species' responses to climate change: movement (shift their geographic ranges to environmental conditions within which they are able to maintain populations), adaptation (in terms of evolutionary change or of physiological acclimatization), or face extinction (Holt, 1990).Ecological niche modeling is an efficient and widely disseminated approach in conservationist studies to predict species' spatial and evolutionary responses to the effects of global climate change via identification of their environmental requirements (Soberón & Nakamuna, 2009).In this paper, we apply an ecological niche model by the MaxEnt algorithm to modeling distributional responses in the face of climate for the two subspecies of Melipona quadrifasciata Lepeletier.
This species is one of the best-known social bee found in Brazil, acting as a pollinator of several plant species and playing an important ecological role in the ecosystem (Ramalho et al., 2004).In nature the species constructs the nest inside trees hollows (Michener, 2007), so it is sensitive to forest fragmentation, especially when nesting sites are destroyed and availability of trophic resources is declined, compromising the maintenance of natural colonies.
M. quadrifasciata has been divided into two subspecies: M. q. quadrifasciata and Melipona q. anthidioides Lepeletier (Moure et al., 2007).The geographic distribution of each subspecies seems to be very distinct.M. q. quadrifasciata occurs in the Southern states of Brazil, from Rio Grande do Sul to the southern São Paulo state, mainly in cold regions, whereas M. q. anthidioides is found in habitats with higher temperatures from northeastern São Paulo state to the northern part of Bahia, and westwards to the western tip of Minas Gerais and central region of the Goiás state (Kerr, 1976).The main morphological difference between the two subspecies consists in the presence of continuous tergal stripes from the 2nd to the 5th segment in M. q. quadrifasciata and interrupted stripes in M. q. anthidioides (Schwarz, 1948).Nonetheless, there is a hybridization zone between the central region of São Paulo state and the southern state of Minas Gerais, where intermediate patterns of tergal stripes can be found (Kerr, 1976).Furthermore, populations with an abdominal pattern identical to that of M. q. quadrifasciata could be found in the northeastern part of Bahia state in warm regions with altitudes ranging from 500 to 700 m (Batalha-Filho et al., 2009).
Recently, molecular studies have been conducted on these species to identify and substantiate the maintenance of division into subspecies.While some studies have provided evidence that it is possible to recognize molecular differences between the two M. quadrifasciata subspecies, confirming the usefulness of the yellow metasomal stripes to identify them (Waldschmidt et al., 2000;Moretto & Arias, 2005;Souza et al., 2008;Tavares et al., 2013), Nascimento et al. (2010) did not reveal genetic structuring of M. quadrifasciata in function of the tergite stripe pattern.Similar results were reported by Batalha-Filho et al. (2009) through PCR-RFLP, and in phylogeographic studies of these bees (Batalha-Filho et al., 2010).
The maintenance of distinct tergal strip patterns between M. q. quadrifasciata and M. q. anthidioides at different geographic regions together with the occurrence of a hybridization zone in the central region suggest that ecological features, such as adaption to distinct climatic conditions may take place (Kerr, 1976).Thus, the possibility of having different responses in geographic range shifts to future climate scenario would be expected.In view of this scenario, in the present study, we investigated the effects of climate changes on the distribution of the two M. quadrifasciata subspecies in Brazil to answer two main questions: (i) Are there differences in geographic shift patterns of the two subspecies to future climate scenario?(ii) Is there an overlap of suitable habitats between them and, if so, how much will it be?Therefore, the present study clarifies questions regarding geographic distribution, and climate niche overlap of M. quadrifasciata subspecies.

Material and Methods
Occurrence data of subspecies of M. quadrifasciata were obtained from the following databases: the Global Biodiversity Facility -GBIF (www.gbif.org)and SpeciesLink (www.splink.cria.org.br).The records that were not previously classified at the subspecies level were identified through available images from the collections and, when this was not possible, were classified according to the distribution map of Batalha-Filho et al. (2009).We excluded data that: (i) came from meliponaries (to avoid over-representation of some areas); (ii) were not georeferenced; and (iii) were impossible to classify to subspecies level.Thus, 137 records were obtained for M. q. anthidioides and 68 records were obtained for M. q. quadrifasciata (Fig 1).
The environmental variables were obtained from the WorldClim database (Fick & Hijmans, 2017) at a resolution of 10 arc-minute grid cells, what corresponds to 0.16 degrees or approximately 18 km in the tropics.The current climate model was based on the time series produced from 1970-2000.Future climate models were developed by the IPCC5-AR5 Fifth Assessment Report (Field et al., 2014) in the scenario of greenhouse gases RCP8.5 to 2100.This scenario is considered the most pessimistic projection, predicting a global temperature rise by about 5 to 6 °C by 2100 and increase of extreme events such as heat waves and very intense rains (Gent et al., 2011).It was chosen because we are currently on a trajectory closer to the RCP8.5 scenario than to the more moderate scenarios (Peacock, 2012).
Considering the climatic conditions of the current subspecies distribution (tropical and subtropical environments and respective transition zones with seasonal temperature and high humidity (Michener, 2007), 11 climatic variables were pre-selected.In order to avoid overfitting and reducing the predictive power of the model, variables that were correlated were excluded using Variance Inflation Factor -VIF (Dormann et al., 2013).We used the function vifstep included in the package usdm (Naimi, 2015), using 10 as threshold.The final dataset of variables selected for the models were: annual average temperature (BIO1), isothermality (BIO3), annual temperature rate (BIO7), precipitation of the wettest month (BIO13), precipitation of the driest month (BIO14), and seasonality of precipitation (BIO15).
A maximal entropy model (MaxEnt) was used to model the current distributions of subspecies of M. quadrifasciata and to design their future habitat suitability.MaxEnt is a machine learning approach which is based on maximum entropy algorithm (Phillips et al., 2006).MaxEnt has been show an adequate and widely method for modelling presence-only data (Elith et al., 2006).This approach uses the presence-only data and background points (pseudoabsences), relating environmental predictors to construct suitability maps ranging from 0 (unsuitable) to 1 (suitable).
In both subspecies we allowed linear, quadratic, product and hinge features between the habitat suitability values and each covariate (Phillips & Dudı´k, 2008).
The projections were repeated 10 times, each selecting a different random sample of 80% when verifying the accuracy of the model in relation to the remaining 20% (Phillips et al., 2006).The MaxEnt algorithm performs better when the study area for the calibration model does not include areas outside the occurrence of species records (Elith et al., 2011).Therefore, considering the current distribution of the two subspecies, the obtained models were limited to the Atlantic Forest, Pampa, and Cerrado biomes.The models' performances were evaluated by Area Under the receiver-operator curve (AUC).To reduce sampling bias, a polarization layer was constructed, consisting of a map that incorporates sampling effort, that is, making the model less important for areas with high density of occurrence points and increasing the importance of regions with similar environmental characteristics, but which were not sampled (Kramer-Schadt et al., 2013).
To characterize potential changes in the distributions of the two subspecies, we calculated the proportions of occupied grid cells in actual and projected scenarios in the Atlantic Forest, Pampa, and Cerrado delimitation.In order to produce the presence/absence maps to calculate the differences in occupied cells we first calculated the threshold which maximized the TSS.We used the functions included in the Presence Absence package (Freeman & Moisen, 2008).
The niche similarities and overlap between the two subspecies in present and future climatic scenarios was verified using the method proposed by Broennimann et al. (2012), in a grid-environment space.This method uses core density functions to calculate the smoothed density of the number of occurrences and the available environments along the environmental axes of the Principal Components Analysis (PCA) and, based on these values, an occupancy index is estimated.Afterwards, the Schoener D-index was calculated, quantifying the niche overlap between the two subspecies (Warren et al., 2008;Broennimann et al., 2012).
Finally, two randomization tests were performed to evaluate the equivalence and niche similarity between the two subspecies.The equivalence test assesses whether niches are indistinguishable (Warren et al., 2008).For this test, pseudo-replications were created, grouping the geographical areas of occurrences of the two subspecies and randomly dividing them into two groups; however, original sample sizes were maintained.For each of the pseudo-replications, we calculated the D-index and then contrasted the original D-value observed with a null distribution of 100 pseudoreplicated D-values.The hypothesis of niche equivalence was rejected if the probability of observed D falling in the null distribution was less than 0.05 (p < 0.05).
The niche similarity test was used to question whether the occupied environmental space in one range is more similar or more different than the occupied environmental space in the other range than would be expected at random.For this test the distribution of the two subspecies in one interval was overlapped with the distribution in the other interval, randomly assigning a new location to each occurrence in the other range and, for each pseudo-repetition, index-D was calculated.This procedure was performed 100 times (from M. q. anthidioides to M. q. quadrifasciata and from M. q. quadrifasciata to M. q. anthidioides) to generate two new null distributions of D-values.The hypothesis of niche similarity was rejected if the probability of observed D falling in the null distribution was less than 0.05 in a two-column test (p < 0.05) (Broennimann et al., 2012).

Results
The AUC values obtained for each model exceeded 0.7, with the highest value (0.785) occurring in the future model of M. q. anthidioides, which indicates a good agreement of the models and good predictive performance.
Our distribution model for the future scenarios of M. q. anthidioides suggested a clear trend of reduction in habitat suitability in the Cerrado biome, and a potential increase of suitable climatic conditions in the Atlantic Forest, especially on the eastern coastal line from the state of Santa Catarina to the state of Espírito Santo.Furthermore, the model predicted a shift from north to south towards the Pampa biome (Fig 2 A, B).Meanwhile, the resulting distribution of M. q. quadrifasciata in future climatic scenarios showed a potential decrease in habitat suitability in almost all of its current distribution, with an increase of suitable landscapes only in the coastal region of Bahia (Fig 2 C, D).In addition, it is predicted to be an increase in suitability of M. q. anthidioides towards the potential distribution area of M. q. quadrifasciata (Fig 2 B, D).
The percentage of appropriate area occupation (with threshold of 0.5 for both subspecies) was 30.35% for M. q. anthidioides in the current model and increased to 31.75% in the future model.For M. q. quadrifasciata, the percentage of suitable area pixels was 16.36% in the current model and decreased to 13.74% in the future model.
The current and future niches of the subspecies in the environmental space are defined by the first two axes The current and future PCAs showed a displacement of the centroid of the M. q. anthidioides niche (area with higher probability density of occurrence in the environmental space) towards BIO3, BIO14, and BIO15 compared to the centroid of the M. q. quadrifasciata niche (Fig 3 A, B  .(E) Niche similarity of M. q. quadrifasciata to M. q. anthidioides, and (F) niche similarity of M. q. anthidioides to M. q. quadrifasciata.

Discussion
According to our potential distribution models, temperature and precipitation (from the driest month and seasonality) are determinant factors for changes in distributions of the two subspecies studied.Several studies on foraging activities of Melipona species showed speciesspecific characteristics of flight in response to the variations in abiotic conditions (e.g.Teixeira & Campos, 2005), allowing interference on the geographic and ecological distribution of each species (Pereboom & Biesmejer, 2003).These studies conclude that air temperature and relative humidity directly affect the flight activity of Melipona species with relatively large body sizes.
Scientific models of climate change for future climatic conditions indicate for the Cerrado biome that, although rainfall tends to increase mainly in the form of more intense and extreme rainfall events, periods of drought will be longer and summers will become warmer with an increase in temperature of up to 4 °C by 2100 (Juras, 2008).According to our models, as a consequence of the temperature increase and drought period predicted for the Cerrado, M. q. anthidioides and simulated niche overlaps (grey bars).(E) Niche similarity of M. q. quadrifasciata to M. q. anthidioides, and (F) niche similarity of M. q. anthidioides to M. q. quadrifasciata.
will undergo a notable reduction in potential habitats in this biome, but will have a potential increase of geographical distribution in the Atlantic Forest, along almost the entire eastern coastline of Brazil, while M. q. quadrifasciata will suffer a reduction of adequate habitats in almost all of its current geographic distribution.This finding reinforces the evidence of the existence of specific responses to variations in climatic conditions, as proposed by Teixeira and Campos (2005).
Furthermore, our results indicate that the two subspecies showed clear differences in geographic shift patterns, supporting the maintenance and acceptance of the division of M. quadrifasciata into two subspecies.
Regarding the disjunct populations with continuous bands similar to that of M. q. quadrifasciata in the northeast of Brazil, our distribution models showed similar patterns in responses to climatic changes with M. q. anthidioides.Waldschmidt et al. (2000) and Batalha-Filho (2010) showed that these northern populations with continuous tergal stripes and colonies of M. q. anthidioides were genetically similar, indicating that the existence of similar tergal band patterns does not determine that this group belongs to the subspecies M. q. quadrifasciata found in the southern part of Brazil.Furthermore, the similar pattern in responses to climate changes detected in our study, added to the low genetic variability among colonies of the two M. quadrifasciata groups found by these authors, suggests that both belong to the same subspecies, M. q. anthidioides.Perhaps, the pattern of continuous bands of disjunct populations would be the result of the more recent differentiation within this group, resulting in a group adapted to higher temperatures and lower precipitation that characterizes the region of the northeast occupied by these populations.
According to our distribution models, M. q. anthidioides will show an increase in distribution suitability towards the Pampa biome.Both subspecies show low current habitat availability in this biome, probably due to lack of adequate nesting sites, since Melipona species build their nests in hollow trunks of trees that are not abundant in this region.This biome is formed by four main groups of natural field vegetation: Campanha Plateau, Central Depression, Sul-Rio-Grandense Plateau, and Coastal Plain, with herbaceous and shrub vegetation predominating and trees restricted to 'capões' and riverbanks (IBGE, 2004).For this reason, the migration of M. q. anthidioides to this biome is unlikely, although in the future this region presents favorable climatic conditions for the subspecies.
Our distribution models clearly indicated that the future climate scenario will be more favorable to M. q. anthidioides when compared to M. q. quadrifasciata, as the first presented greater potential to expand its distribution through migration.This migration, according to our results, would lead to a 0.7% increase in the overlap of the potential distribution between the two subspecies.The increase in overlap, in turn, may lead to an increase in the size of the hybridization area currently existing between the two subspecies (Batalha-Filho et al., 2009).
Hybridization between populations may allow alleles from one genetic background to integrate into another if favored by selection (Riesenberg, 2003).However, the net outcome of inter-population hybridization can be affected by several factors, such as the level of divergence between the populations and the rate of inbreeding in the populations (Whitlock, 2000).In the present hybridization zones of the two subspecies of M. quadrifasciata, there are individuals with dorsal bands similar to one subspecies but with genetic markers of the other subspecies (Waldschmidt et al., 2000;Batalha-Filho et al., 2009).In view of the future scenarios of the potential geographical distribution of the two subspecies (increase of suitable habitats of M. q. anthidioides, potential decrease in distribution for M. q. quadrifasciata, and increased climate niche overlap) it is possible that M. q. quadrifasciata, which is already in the red list of the endangered fauna of the state of Rio Grande do Sul (Fontana et al., 2003), will be at high risk of loss.
Due to the variation of some genetic markers, as RAPD and RFLP in the cytochrome b gene, social bees tend to have a differentiated genetic predisposition between their colonies in relation to their foraging behavior for example (Oldroyd et al., 1992).The loss of M. q. quadrifasciata may lead to loss of differential responses presented by this subspecies, and, consequently, result in a reduction in the genetic variability of M. quadrifasciata groups.Araújo et al. (2000) showed that populations of Melipona are particularly sensitive to high genetic drift due to homozygosity in X 0 sex determination locus, which generates diploid males that are either unviable or sterile, putting at risk the maintenance of the species.
In this sense, to avoid that M. q. quadrifasciata goes extinct, conservation actions aiming to increase suitable areas for the expansion of this subspecies are urgently needed.

Fig 1 .
Fig 1. Spatial distribution of records of the two Melipona quadrifasciata subspecies in Brazil.Grey circles represent the records of M. q. anthidioides, and black squares the M. q. quadrifasciata records.
of the PCA, which explained 76.7% of the original climatic variation.The subspecies showed little overlap of the current niche (Schoener's D = 0.399, Fig 3 D), and an increase in future niche overlap (Schoener's D = 0.406, Fig 4 D).

Fig 3 .
Fig 3. Current climatic niches occupied by the two subspecies of Melipona quadrifasciata.Niche occupied by M. q. anthidioides (A) and M. q. quadrifasciata (B) along the two-first axes of the PCA (see (C) for details).Grey shading represents the density of the occurrences of the subspecies by cell.The solid and dashed contour lines illustrate, respectively, 100% and 50% of the available environment.(C) Contribution of the climate variables to the first-two axes of the PCA (bio1: annual mean temperature, bio3: isothermality, bio7: annual temperature rate, bio13: precipitation of the wettest month, bio14: precipitation of the driest month, and bio15: precipitation seasonality), and contribution of the two-first axes to data variation.(D) Histogram of the observed niche overlap D (D = 0.399) (bars with a diamond) and simulated niche overlaps (grey bars).(E) Niche similarity of M. q. quadrifasciata to M. q. anthidioides, and (F) niche similarity of M. q. anthidioides to M. q. quadrifasciata.

Fig 4 .
Fig 4. Future climatic niches occupied by the two subspecies of Melipona quadrifasciata.Niche occupied by M. q. anthidioides (A) and M. q. quadrifasciata (B) along the two-first axes of the PCA (see (C) for details).Grey shading represents the density of the occurrences of the subspecies by cell.The solid and dashed contour lines illustrate, respectively, 100% and 50% of the available (background) environment.(C) Contribution of the climate variables to the first-two axes of the PCA (bio1: annual mean temperature, bio3: isothermality, bio7: annual temperature rate, bio13: precipitation of the wettest month, bio14: precipitation of the driest month, and bio15: precipitation seasonality), and contribution of the two-first axes to data variation.(D) Histogram of the observed niche overlap D (D = 0.406) (bars with a diamond)and simulated niche overlaps (grey bars).(E) Niche similarity of M. q. quadrifasciata to M. q. anthidioides, and (F) niche similarity of M. q. anthidioides to M. q. quadrifasciata.