Evidence for morphological and genetic structuring of Plebeia flavocincta (Apidae: Meliponini) populations in Northeast Brazil

Maia, Ulysses Madureira; Santos Júnior, José Eustáquio dos; Molina, Michele; Galaschi-Teixeira, Juliana Stephanie; Carvalho, Airton Torres; Miranda, Leonardo de Sousa; Imperatriz-Fonseca, Vera Lucia; Oliveira, Guilherme; Giannini, Tereza Cristina

doi:10.3389/fevo.2022.1057624

ORIGINAL RESEARCH article

Front. Ecol. Evol., 25 November 2022

Sec. Evolutionary and Population Genetics

Volume 10 - 2022 | https://doi.org/10.3389/fevo.2022.1057624

Evidence for morphological and genetic structuring of Plebeia flavocincta (Apidae: Meliponini) populations in Northeast Brazil

$\r\nUlysses Madureira Maia,*&#x;$ Ulysses Madureira Maia^1,2*†

José Eustáquio dos Santos Júnior^3†

Michele Molina^2,4†

Juliana Stephanie Galaschi-Teixeira^3†

Airton Torres Carvalho^5†

Leonardo de Sousa Miranda^2,6†

Vera Lucia Imperatriz-Fonseca^7†

Guilherme Oliveira^2,4†

Tereza Cristina Giannini^1,2†

¹Instituto de Ciências Biológicas, Universidade Federal do Pará, Belém, Brazil
²Instituto Tecnológico Vale, Belém, Brazil
³Instituto de Genética e Bologia Molecular, Universidade Federal do Pará, Belém, Brazil
⁴Departamento de Genética, Ecologia e Evolução, Instituto de Ciências Biológicas, Universidade Federal de Minas Gerais, Belo Horizonte, Brazil
⁵Departamento de Ciência Animal, Universidade Federal Rural do Semi-Árido, Mossoró, Brazil
⁶Lancaster Environment Centre, Lancaster University, Lancaster, United Kingdom
⁷Instituto de Biociências, Universidade de São Paulo, São Paulo, Brazil

Tropical dry forests are increasingly threatened by human activities. In Northeast Brazil, the Caatinga, an area of tropical dry forests surrounded by wetlands, is considered globally unique among these habitats. In this region, the stingless bee Plebeia flavocincta is found in a variety of environmental, ecological, and demographic conditions. We aimed to characterize P. flavocincta populations within its natural range through wing geometric morphometrics and mitochondrial DNA (mtDNA) analyses. The characterization of population variability can clarify whether the species is morphologically and genetically diverse and whether populations are morphologically and genetically structured. We analyzed 673 samples by wing morphometry and 75 by cytochrome-b assays. Our results revealed P. flavocincta is genetically and morphologically diverse and populations are morphologically and genetically structured. Despite the differentiation between the two most morphologically distant populations, we verified a large overlap of morphological variation between all populations. The genetic analysis showed that the haplotypes were geographically structured into six clusters, four of which were located in coastal areas, and the remaining two in the drier inland region. The characterization of P. flavocincta populations is an important step toward decision-making in programs for the protection, management, and sustainable use of this stingless bee in local breeding efforts.

Introduction

Arid and semiarid regions cover more than 30% of Earth’s land surface (Okin et al., 2006). Within these regions, tropical dry forests are estimated to cover more than one million square kilometers, with more than half in South America (54%), but most are under threat (Miles et al., 2006). In Northeast Brazil, tropical dry forests, locally known as the “Caatinga,” are globally unique, presenting high levels of endemism (Silva et al., 2017; Fernandes et al., 2020). In this region, the local biota is subjected to long periods of high heat and irregular water availability for several months each year. In recent decades, this region of Brazil has undergone a rapid and extensive environmental change caused by the over exploitation of its natural resources (Silva et al., 2017), which are naturally limited due to the region’s climatic conditions (Krol and Bronstert, 2007). Furthermore, future climate projections show large temperature increases and rainfall reductions, suggesting the occurrence of more frequent and intense droughts and a tendency toward desertification in this region (Marengo et al., 2017). Studies aiming to analyze the genetic and phenotypic variability of populations are urgently needed as they can assess the degree of inbreeding, gene flow, and genetic diversity among and within populations (Keyghobadi et al., 2005; Wagner, 2012; Carvalho et al., 2019). Given the increasing loss and modification of natural habitats, in the long term, loss of genetic variation is expected to affect the ability of populations to respond to changing selection pressures (Cheptou et al., 2017).

Changes in natural habitats and species composition can impact ecosystem services and compromise their long-term continuity (González-Varo et al., 2013). One important ecosystem service is pollination, as this service facilitates reproduction in the vast majority of flowering plant species (Potts et al., 2016). Pollination success is closely related to the presence and availability of pollinators (Campbell et al., 2018). Most species of tropical flowering plants depend on insects and other animals as pollinators (Ollerton et al., 2011). Among pollinating animals, bees are considered the most important (Freitas et al., 2009; Ollerton, 2017). Many bee species native to the seasonally dry tropical forests of northeast Brazil, notably species of stingless bees (Meliponini tribe) (Michener, 2013), are exploited for honey production by local beekeepers (Jaffé et al., 2015).

Plebeia flavocincta (Cockerell 1912) is a small (∼4.5 mm) stingless bee with a distribution that is restricted to Northeast Brazil. Its distribution covers both coastal areas and seasonally dry tropical forests (Maia et al., 2020). In this region, the species occurs in a variety of environmental, ecological, and climatic conditions. P. flavocincta can nest in exotic trees and urban areas (Ribeiro and Taura, 2019) and is also considered a docile species, which facilitates its rearing in hives for honey, geopropolis, and pollen (Maia et al., 2015). In addition, it can be used as a pollinator in greenhouses to increase crop production (Imperatriz-Fonseca et al., 2017).

Several studies on pollinators such as bees have demonstrated that tools such as geometric morphometrics and mitochondrial DNA (mtDNA) analysis are robust for detecting subtle variations in populations (Bonatti et al., 2014; Rattanawannee et al., 2017). The geometric morphometry of the wings is useful for detecting small morphological variations that are imperceptible by traditional morphometry (Francoy et al., 2011; Nogueira et al., 2019). The CytB mtDNA gene is useful not only in identifying cryptic bee species but also in assessing intraspecific distribution patterns (Françoso and Arias, 2013; Santos-Júnior et al., 2015, 2019). The application of complementary methods to analyze population variability has shown promising results in understanding the distribution of species and in resolving taxonomic problems (Bonatti et al., 2014; Francoy et al., 2016; Galaschi-Teixeira et al., 2018).

Considering that organisms living in regions with high ecological variation tend to present corresponding levels of biological variation, we aimed to characterize P. flavocincta populations within their natural range through wing geometric morphometry and mtDNA analysis. The characterization of population variability can clarify whether the species is morphologically and genetically diverse and whether populations are morphologically and genetically structured.

Materials and methods

Study area

Within its extension, Northeast Brazil is divided into four biomes (Amazon, Cerrado, Caatinga, and Atlantic Forest). The annual accumulated precipitation does not exceed 700 mm in some areas of the seasonally dry tropical forest, which is frequently affected by long periods of drought (Alvares et al., 2013). In contrast, there are areas in the coastal region where the annual rainfall is more than 1500 mm, and the climate is predominantly a tropical humid climate (Awange et al., 2016; Oliveira et al., 2017).

We collected samples in coastal areas and areas of seasonally dry tropical forest (Caatinga) (Figure 1). The coastal areas contain a wide variety of vegetation native to the coastal region, including deciduous seasonal forest, semideciduous seasonal forest, dense ombrophilous forest, open ombrophilous forest, mixed ombrophilous forest, and associated ecosystems such as sandbanks (restingas), mangroves, and altitude fields (Oliveira-Souza et al., 2018). The Caatinga areas include species typical of shrub communities, with vegetation approximately 8 m in height and predominantly consisting of sparsely distributed small trees and shrubs, in addition to herbaceous patches that develop only during the rainy season (Tavares-Dasmasceno et al., 2017).

FIGURE 1

Figure 1. Plebeia flavocincta sampling locations across Northeast Brazil over a map with the two main biomes of the region (Caatinga and Atlantic Forest) showing the two methods used (wing morphometry and molecular analysis). Only one type of analysis was used in some locations.

We collected workers of P. flavocincta from natural nests in 18 locations within their natural range. For the study of the geometric morphometry of the wings, due to differences in the number of nests studied between the localities, we grouped some nearby localities (∼100 km) into microregions (Table 1) according to the Instituto Brasileiro de Geografia e Estatística (IBGE) (Instituto Brasileiro de Geografia e Estatística, 2017). This also reduced the quantity of localities and increased the number of samples per region, increasing the robustness of the analyses.

TABLE 1

Table 1. Quantity of studied nests of Plebeia flavocincta (Apidae: Meliponini) at 18 locations (15 nests wing geometric morphometry analysis and 16 nests for mitochondrial DNA analysis).

Morphometric analysis of wings

We collected at least 10 workers of P. flavocincta from 71 nests. To analyze the venation pattern of the wings, the right forewing was removed with the aid of tweezers and later mounted between the blade and a coverslip with a brush. Images of these wings were captured by a digital camera attached to a stereomicroscope (Olympus) and stored on a computer. Subsequently, a “.tps” file was created with the images of the wings using tpsUtil software version 1.40 (Rohlf, 2008b), and 13 anatomical landmarks were manually marked at the vein junctions of the cells of each wing using tpsDig software version 2.12 (Rohlf, 2008a) (Figure 2). The images were aligned by Procrustes (Bookstein, 1992) and the average configuration of the wings of bees from the same nest were used as a comparative parameter. We performed principal component analysis (PCA) as an exploratory analysis to characterize shape variation using the relative Cartesian coordinates of each landmark after alignment. Next, we performed canonical variate analysis (CVA) to determine the characteristic structure that best distinguished the groups. We analyzed groups by assigning colonies to geographic origins (Table 1). These analyses were performed using the program MorphoJ version 2.0 (Klingenberg, 2011). To verify if the most morphologically similar populations were also the closest geographically, we also calculated the Mahalanobis distances for all pairs using the pooled within-group covariance matrix for all the groups jointly and correlated them with the geographical distances between the populations using a Mantel test. To graphically visualize the morphological proximities between populations, we created a morphological proximity dendrogram using MEGA version 5 (Tamura et al., 2011).

FIGURE 2

Figure 2. Plebeia flavocincta (A) and its wing (B) showing the 13 anatomical landmarks used (Photos: UMM).

Molecular analyses

For mtDNA analysis, we used P. flavocincta individuals from 75 nests (one individual per nest) from 16 locations within the natural range (Table 1). The laboratory procedures described below were carried out at the Instituto Tecnológico Vale Desenvolvimento Sustentável, Belém, Brazil.

Deoxyribonucleic acid extraction, amplification, and sequencing

Genomic DNA was extracted from one worker per colony using the DNeasy Blood & Tissue Kit (Qiagen) following the manufacturer’s instructions. The Cytochrome B (CytB) region of mtDNA was amplified using the MtD26/MtD28 primer pair (Simon et al., 1994). PCR was performed with a final reaction volume of 25 μl comprising 5–7 μl of template DNA, 5 μl of PCR buffer (5 ×), 2 μl of dNTPs (2 mM), 2.4 μl of MgCl2 (25 mM), 0.5 μl of each primer (10 mM), 2 μl of TBT (trehalose, bovine serum albumin, and polysorbate-20), and 0.2 μl of Go Taq Flex Polymerase (Promega 5 U/μl). Ultrapure water q.s.p. was used to reach a final volume of 25 μl. All PCR amplifications were performed according to Françoso et al. (2016) on a Mastercycler Pro. For amplification with the primers described by Simon et al. (1994), we used the same PCR cycling conditions for the barcode region of DNA. The PCR products were run on a 1.2% agarose gel and visualized with a UV light transilluminator. Subsequently, the amplicons were purified as described in Vasconcelos et al. (2021), and the sequencing reactions were performed with the same primers (forward and reverse) used in the amplification and the amplicons and the reagents of the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) according to the manufacturer’s instructions. Sequencing was performed with an ABI3730 DNA Analyzer (Applied Biosystem).

The sequences obtained were analyzed using the SeqA 6 program (Applied Biosystems) and later edited, checked, and aligned using Geneious Prime 2019.2.3 software (Biomatters). For alignment, MUSCLE (Edgar, 2004) was used. DnaSP v.5.10 software (Librado and Rozas, 2009) was used to separate the haplotypes into their respective groups and determine the number of haplotypes (h) and haplotypic diversity (Hd). The haplotype network was built using NETWORK v. 4.6 (Polzin and Daneshmand, 2003), using the median joining algorithm (Bandelt et al., 1999). The mean number of inter- and intrapopulation variable sites was calculated using the MEGA program version 5 (Tamura et al., 2011). To verify significant correlations between the genetic distances and geographic distances of the studied locations, a Mantel test was performed with the aid of TFPGA software (Miller, 1997).

Population differentiation was detected using Bayesian clustering algorithms in the program GENELAND ver. 4.0.0 (Guillot et al., 2005). First, we performed two runs of 1,000,000 Markov chain Monte Carlo (MCMC) iterations with a thinning interval of 100 and a maximum number of populations of 10 (K = 10). We used correlated allele frequencies and geographic coordinates to parameterize all runs. We inferred K as the modal number of genetic groups estimated among the best of the 1,000,000 iterations for the two runs. Then, we performed two independent runs with K fixed to the number of populations previously inferred.

Demographic inference

A Bayesian Skyline Plot (BSP) reconstruction was implemented to evaluate population size dynamics over time using the BEAST 2.6.6 program (Bouckaert et al., 2014). A CytB substitution rate of 1.9% per lineage per million years per generation (1 year) was adopted [see (Miranda et al., 2017)]. The evolutionary model was inferred by BEAST Model Test (all reversible). The following parameters were employed: a strict clock, 20 million generations, parameter sampling every 1,000 generations of the MCMC analysis, and a 10% burn-in. Convergence between runs and the performance of the analysis were checked using Tracer 1.7.1 and the plots were constructed with the aid of the same program (Rambaut et al., 2018).

Results

Morphological characterization of populations through geometric morphometry

Our data revealed at least 22 measurements of relative deformations in the geometric morphometry of the wings, according to the Procrustes residuals generated for PCA (Figure 3). Among such measures, the first 14 explained 91.46% of the variability in the entire data set, with PC1 and PC2 explaining 22.15 and 13.92% of the variation in the sample, respectively. In the scatter plot of PCA results (Figure 3), the horizontal axis (PC1) shows differentiation with large overlapped variation space among these groups. Axis 2 (PC2) showed that there are no differences between the pairs of populations since they completely overlapped.

FIGURE 3

Figure 3. Principal component analysis of the 13 anatomical landmarks marked at the wing rib junctions of Plebeia flavocincta.

Canonical variate analysis using the covariance matrix of all groups showed eight axes that maximized the separation of these populations (Table 2). Among the eight axes, the first two axes (CV1 and CV2) explained 55.71% of the variation in the data set. The characteristic distribution of points that best showed the separation of P. flavocincta groups can be seen in the scatter plot (Figure 4). Although separation was found between the RN coastal, Ibiapaba, and Maciço, there was an overlap of morphological variation of the other pairs between them.

TABLE 2

Table 2. Variation of the axes of canonical variables.

FIGURE 4

Figure 4. Canonical variate analysis of the 13 anatomical landmarks marked at the wing rib junctions of Plebeia flavocincta.

Considering the squared Mahalanobis distances (Table 3) and Procrustes distances between populations (Table 4), the greatest morphological proximity (shortest distance) occurred between the Cariri and Potengi populations, while the greatest distance was observed between the RN coastal and PE coastal populations. The RN coastal population and, primarily, PE coastal population showed the greatest distances from the populations of all other locations (Table 3). The Mantel test based on the Mahalanobis distances and geographic distances showed no correlation between morphological and geographic distances (r = 0.060, p = 0.321).

TABLE 3

Table 3. Squared Mahalanobis distances between the centroids of the analyzed regions.

TABLE 4

Table 4. Procrustes distances between regions obtained from the anatomical landmarks of the wings of Plebeia flavocincta.

Based on the squared Mahalanobis distances, a dendrogram of morphological proximity between the studied groups was constructed (Figure 5), in which the coastal populations were morphologically farther apart and formed groups that were more external and isolated from each other than populations from inland areas.

FIGURE 5

Figure 5. Morphological proximity dendrogram based on the neighbor-joining method calculated from the Mahalanobis distances between the groups (regions) studied.

Genetic characterization of populations through mitochondrial DNA

Eighteen haplotypes were found (Figure 6 and Table 5), with Hap-15 being the most frequent (23 colonies). Haplotypes Hap-01, Hap-02, Hap-04, Hap-06, Hap-10, Hap-12, Hap-14, and Hap-18 were less frequent, occurring in only one colony each. Hap-3 was shared by four locations. The haplotype diversity index (Hd) was 0.85. Some localities had a single haplotype that was not shared with any other locality. The clustering approach implemented in Geneland consistently detected six genetic clusters (subpopulations) (Figure 7). These results suggest separation of the four populations along the coast and two other populations in inland areas of Northeast Brazil.

FIGURE 6

Figure 6. Haplotype network showing shared haplotypes of the observed populations obtained using the Median Joining algorithm.

TABLE 5

Table 5. Locations where each haplotype of Plebeia flavocincta was originally found.

FIGURE 7

Figure 7. Map of the locations where P. flavocincta samples were collected in Northeast Brazil for genetic analysis. The colors indicate the six different populations identified by Geneland software (K = 6).

Population size dynamics over time

The BSP did not reveal any substantial fluctuations in the effective population size (Ne) in P. flavocincta in the last 200,000 years (Figure 8).

FIGURE 8

Figure 8. Coalescent Bayesian Skyline Plot (BSP) used to infer demographic history of Plebeia flavocincta population. Black horizontal line shows median BSP estimate and gray area shows upper and lower 95% highest posterior density limits.

Discussion

Both the morphometric and molecular analyses were useful for inferring structure among the studied populations. Although there was substantial overlap between methods, results from the morphometric analysis showed that the most geographically distant population, on the coast of PE coastal, corresponded to the most differentiated population. Despite the discrimination between the two most geographically distant populations (in this case, PE coastal and Ibiapaba), we could verify the occurrence of an overlap of morphological variation between them, formed by the other populations. Moreover, the genetic analysis showed that the haplotypes were geographically structured in six clusters, with four of them located in coastal areas and two others in inland areas.

Geographically closer populations, such as those on the coast, were more morphologically distinct from each other. On the other hand, some more geographically distant populations within the Caatinga areas were morphologically similar. An explanation for this finding may be the large extent of the northeast, with areas of large thermal fluctuations, variability in vegetation cover, with deforested areas and mountain ranges acting as dispersion barriers (Jaffé et al., 2019). Isolated mountain ranges and high plateaus (up to 1000 m) are spread across the Caatinga areas. On the borders of the State of Ceará, for example, there are mountain ranges that can prevent dispersion and promote local adaptation of populations. In addition, coastal areas have variable amounts of forest cover that extends over a wide latitudinal and longitudinal strip, which has favored high levels of diversity and endemism of both the fauna and flora (Ribeiro et al., 2009). A latitudinal effect on the distribution of genetic variation has already been observed for Melipona subnitida Ducke, 1910, a species of stingless bee that has an area of occurrence similar to that of P. flavocincta (Jaffé et al., 2019). Still, the isolation of populations in the coastal region of the Northeast may be a reflection of the fragmentation of the Atlantic Forest. Another plausible explanation is that gene flow was facilitated in Caatinga areas, probably by the transport of colonies originating from meliponiculture practices (Jaffé et al., 2016). Many species of stingless bees are commercially reared for multiple purposes, such as honey production, crop pollination and recreation. However, due to this process, hives can be moved for great distances (Santos et al., 2022). In this sense, the differentiation of the coastal population of the Northeast may be related to the fact that these bees are less intensively managed than those in the Caatinga areas, since beekeeping practices can play an important role in the formation of genetic differentiation between bee species (Jaffé et al., 2016).

The number of haplotypes and haplotype diversity index (Hd) obtained for P. flavocincta were higher than those observed for other stingless bees, such as Mourella caerulea (Friese 1900) (Hd = 0.74) with six haplotypes (Galaschi-Teixeira et al., 2018); M. subnitida (Hd = 0.79) 11 haplotypes (Bonatti et al., 2014); Partamona helleri (Friese, 1900) (Hd = 0.55) with 10 haplotypes (Brito and Arias, 2010); and P. remota (Holmberg, 1903) (Hd = 0.54) with 15 haplotypes (Francisco and Arias, 2010). The presence of unique mitochondrial haplotypes may be associated with the process of colonization of new areas by one or a few colonies of stingless bees (Miranda et al., 2016). This pattern of unique haplotypes has already been observed in M. subnitida (Bonatti et al., 2014; Barbosa et al., 2019). Moreover, a study using thousands of genetic markers of M. subnitida showed two very distinct populations at the extremes of its distribution area and two populations with greater variability in the central region, suggesting that the species colonized the peripheral regions more recently, while the populations from the central region had more time to develop greater genetic variability (Jaffé et al., 2019). The low dispersal capacity of daughter colonies during the swarming process due to dependence on the mother colony, as well as the queen’s philopatry, may help explain these patterns in stingless bees (Engels and Imperatriz-Fonseca, 1990; Nogueira-Neto, 1997). Other studies using population genetic approaches with different species of stingless bees also found similar results (Tavares et al., 2007; Francisco et al., 2008; Batalha-Filho et al., 2010; Brito et al., 2013; Galaschi-Teixeira et al., 2018). Thus, the presence of common haplotypes of P. flavocincta in different and distant locations cannot be explained by current gene flow via females, and the exclusivity of several observed haplotypes can be considered an indicator of population structure. For P. remota, the hypothesis of long-distance migration based on queen philopatry was discarded in favor of human-mediated “colony migration” as a leisure activity (Francisco and Arias, 2010).

Pleistocene climatic instability drove the historical distribution of forest islands in the northeastern Brazilian Atlantic Forest (Silveira et al., 2019). In the wetter areas on the northeast coast, P. flavocincta was likely isolated in the remnants of the Atlantic Forest for long periods, which limited gene flow and allowed the independent accumulation of base substitutions in mtDNA, which may have led to the current differentiation of the populations in these areas. Additionally, according to paleodistribution and palynological validation models, the potential historical distribution of the Caatinga areas was much smaller during the driest and coldest period of the Maximum Glacial, with subsequent expansion with small fluctuations until reaching the current area (Werneck et al., 2011). These fluctuations in suitable habitats may have promoted some isolation among populations. Similarly, models have suggested that P. flavocincta followed these distributional dynamics in the past, and an increase in climatically suitable areas on the northeast coast is expected in the future (Maia et al., 2020).

Our results also revealed patches of other populations within the area attributed to the coastal population of Paraíba and Pernambuco states, a result very similar to that found in areas of climatic refuge (Maia et al., 2020). These areas can be considered the ancestral areas of the populations located in the northern portion of Northeast Brazil (states of Ceará and Rio Grande do Norte). Therefore, this dynamic of expansion of vegetation cover and areas of occurrence seems to explain the current morphological and genetic patterns of the species. However, the Pleistocene climatic changes do not seem to have had a substantial influence on the demographic patterns of the P. flavocincta, since the BSP did not reveal any substantial fluctuations in the effective size of the populations. Although there has been a fluctuation in the climatically suitable areas, this does not have drastically affected the population dynamics. Based on the already demonstrated climatic niche distribution models (Maia et al., 2020), the populations of P. flavocincta probably remained stable even with the decrease in their climatically stable areas. This finding agrees with other studies with stingless bees in the same region (Miranda et al., 2016, 2017).

Here we do not disregard the existence of sampling bias caused by the difficulty in searching for natural nests. In addition to sampling bias, the molecular marker used here only tells the matrilineal history (CytB). However, the present work is an important contribution considering the few studies of population variability with stingless bees in tropical dry forests and aims to help pave the way for a deeper understanding of the evolutionary processes that occur in semi-arid areas. Given its importance as a visitor to native flora and the fact that it is an endemic bee of a threatened dry forest, the present findings may be useful for the development of conservation strategies aimed at P. flavocincta populations. Future studies should prioritize new collections in locations not yet sampled and test how the landscape influences genetic diversity. They must also deepen their ecological knowledge of the species to design conservation programs. Studies involving conservation genetics to maintain the genetic viability of natural populations of the species are also urgently needed.

Conclusion

This is the first study to demonstrate the variability of Plebeia flavocincta populations through the geographic range of the species. We demonstrated that this small pollinator species, occurring in varied environments such as areas along humid coasts and in seasonally dry tropical forests, has morphological and genetic variation to the point of causing population structuring. This likely reflects local adaptation and could be exploited to identify populations more accurately. Furthermore, the observed exclusivity of several haplotypes can be considered an indicator of the gain of new areas and population restructuring. Further work could investigate the ecological or behavioral implications of these findings for this stingless bee species, especially considering the risk of deforestation and climate change in the Caatinga.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/genbank/, OP066721-OP066794.

Ethics statement

Ethical review and approval was not required in accordance with the local legislation and institutional requirements.

Author contributions

UM, VI-F, and TG designed the study. UM and AC collected the specimens. UM, JS-J, MM, and LM performed the analyses. UM, JG-T, GO, VI-F, and TG coordinated the data analysis. UM and JS-J wrote the manuscript. All authors reviewed and edited the manuscript.

Funding

This work was supported by Fundação de Amparo à Pesquisa do Estado do Pará–FAPESPA (doctoral grant, Public Notice 005/2016), the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) scholarship (VI-F, 312250/2018-5), and Instituto Tecnológico Vale (ITV) for infrastructure (project RBRS000603.81).

Acknowledgments

We thank the Universidade Federal Rural do Semi-Árido (UFERSA) and Universidade de São Paulo for logistical support during the study. We also thank Cristina Arias for allowing the internship (UM) to learn genetic procedures, thank Charles Santos, Felipe Contrera, Jeferson Coutinho, Marcia Ribeiro, Orlando Tobias, Tiago Francoy, and William Sabino for suggestions and comments, and also Alistair Campbell for reading and reviewing the text.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alvares, C. A., Stape, J. L., Sentelhas, P. C., Moraes Gonçalves, J. L., and Sparovek, G. (2013). Köppen’s climate classification map for Brazil. Meteorol. Zeitschrift. 22, 711–728. doi: 10.1127/0941-2948/2013/0507

PubMed Abstract | CrossRef Full Text | Google Scholar

Awange, J. L., Mpelasoka, F., and Goncalves, R. M. (2016). When every drop counts: analysis of droughts in Brazil for the 1901-2013 period. Sci. Total Environ. 566, 1472–1488. doi: 10.1016/j.scitotenv.2016.06.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Bandelt, H. J., Forster, P., and Röhl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036

PubMed Abstract | CrossRef Full Text | Google Scholar

Barbosa, M. M., Bonatti, V., Teixeira, J. S. G., Rêgo, M. M. C., and Francoy, T. M. (2019). Genetic characterization of Melipona subnitida stingless bee in Brazilian northeast. Sociobiology 66, 597–601. doi: 10.13102/sociobiology.v66i4.3408

PubMed Abstract | CrossRef Full Text | Google Scholar

Batalha-Filho, H., Waldschmidt, A. M., Campos, L. A. O., Tavares, M. G., and Fernandes-Salomão, T. M. (2010). Phylogeography and historical demography of the Neotropical stingless bee Melipona quadrifasciata (Hym., Apidae): incongruence between morphology and mitochondrial DNA. Apidologie 41, 534–547. doi: 10.1051/apido/2010001