Genetic diversity and population structure of two Euglossini bee species in a host-parasite relationship

Abstract

In the current study, two euglossine species, Exaerete smaragdina and Eulaema nigrita, a cleptoparasite bee and its host, respectively, were used as models to: (i) access the genetic diversity and population structure of both species, sampled along a wide latitudinal range of Atlantic Forest, where the distribution of El. nigrita and Ex. smaragdina co-occurs; (ii) investigate the evolutionary history of these species through the Atlantic Forest, and in a wider scenario, to examine the evolutionary history of these species across others forest domains. Analyses involved males of El. nigrita and Ex. smaragdina sampled through Brazilian territory, including 19 sites in the Atlantic Forest. Bayesian Skyline Plot (BSP) was used to infer possible climate oscillations on population of both species over time. The BSP revealed stability in effective population size for both species in most of the Plio-Pleistocene period. However, BSP results aligned to the starlike configuration in the haplotype network, neutrality test, and population diversity patterns indicated population expansion of the two species during the late Pleistocene. Our findings suggest areas of potential refugia to the climatic oscillations of the Pleistocene in the Atlantic Forest in the Brazilian states of Espírito Santo for El. nigrita and Pernambuco for Ex. smaragdina.

Key words
ecological interaction; evolutionary history; orchid bees; parasitism; Pleistocene

INTRODUCTION

Orchid bees (Hymenoptera: Apidae: Euglossini) are a group of Neotropical pollinators known for their affinity with humid tropical forests (Dressler 1982), which have their origin attributed to Amazon basin (Ramírez et al. 2010). The tribe encompasses five living genera: Eufriesea Cockerell, Euglossa Latreille, Eulaema Lepeletier, Aglae Lepeletier & Serville and Exaerete Hoffmannsegg, (Dressler 1982, Roubik & Hanson 2004). While the first three show free living habits, Aglae and Exaerete are cleptoparasites of other euglossine genera. Specifically, Aglae, a monotypic genus, is an exclusive nest cleptoparasite of Eulaema, and Exaerete species are cleptoparasites of both Eulaema and Eufriesea (Dressler 1982, Roubik & Hanson 2004).

In the cleptoparasitic relationship, the female invades the nest of other species and lays her eggs in the cells of host bees (Garófalo & Rozen 2001, Danforth et al. 2019). Among orchid bees, one of the best known cleptoparasitic relationships is between the species Exaerete smaragdina (Guérin) and Eulaema nigrita Lepeletier. In the first review on orchid bee nesting biology, Zucchi et al. (1969) highlighted that 50% of the cells in a nest of El. nigrita were parasitized by Ex. smaragdina, revealing a high ability of females of this latter species to lay their eggs in the cells of the host species. Subsequently, the cleptoparasitic behavior of Ex. smaragdina in the nest of El. nigrita was described in detail by Garófalo & Rozen (2001).

Eulaema nigrita and Ex. smaragdina have a broad geographic range, occurring from Mexico to southern Brazil (Nemésio 2009). Across the Brazilian territory, both species are also widely distributed and commonly found together in inventories conducted in different vegetation formations (Storck-Tonon et al. 2009, 2013, Silveira et al. 2015, Martins et al. 2018).

It has been suggested that favorable rates of parasitism are usually density-dependent and, consequently, will be enhanced where host populations are large (Wcislo 1987). Despite this, studies revealed that abundances of El. nigrita and Ex. smaragdina are not necessarily positively associated along their distribution. For instance, along the Atlantic Forest domain, it was demonstrated that while the frequency in number of males of El. nigrita, surveyed in inventories, is high towards the South whereas no correlation was found between frequencies of Ex. smaragdina and latitude (Nemésio & Silveira 2006). In fact, in several studies carried out through Atlantic Forest, the relative frequency of Ex. smaragdina varied consistently, in most studies independently from the frequencies of its host (Rebêlo & Garófalo 1997, Aguiar & Gaglianone 2012, Cordeiro et al. 2013). Despite this, a significant positive association between host–parasite ratio, involving both species, and latitude was found by Nemésio & Silveira (2006), who demonstrated that the frequencies of the parasite relative to its host are higher near the equator. They also suggested that there is some trend in frequencies of Ex. smaragdina decrease going further south across the Atlantic Forest (Nemésio & Silveira 2006). While it has been suggested that variations in abundance of some species of orchid bees, probably, reflects variables on smaller spatial scales (Lopes et al. 2022), it is still necessary further investigations on this theme before any generalization.

In this context, studies on the genetic diversity and population structure of these species and their evolutionary history could be helpful to better understand the current abundance and distribution of Ex. smaragdina and El. nigrita across the Atlantic Forest. Furthermore, considering that the relationship between cleptoparasitic bees and their hosts are usually specialized, and the parasitism successful is dependent on both presence and abundance of the hosts, population genetic studies involving parasites and their hosts can be valuable for supporting future conservation and management measures.

Many studies have shown the effects of Pleistocene, a period of climatic changes and geomorphological alterations in the Neotropical region, on the population structure and demography of different organisms in the Atlantic Forest, including different groups of bees (Batalha-Filho et al. 2010, Frantine-Silva et al. 2017). The current literature indicates that paleoclimatic instability, occurring during the Quaternary Period, impacted the Atlantic Forest, shaping the genetic structure of some orchid bee species (López-Uribe et al. 2014, Frantine-Silva et al. 2017). Moreover, it was suggested that species showing narrower physiological tolerance probably experienced less suitable habitats during the Quaternary climatic oscillations (López-Uribe et al. 2014). Thus, orchid bees are excellent models to investigate the effect of climatic oscillations of the Plio-Pleistocene (between 5.3 million years ago–mya and 11,600 years ago) (López-Uribe et al. 2014).

Taking the above into consideration, the aims of the present study were: (i) to investigate genetic diversity and population structure of El. nigrita and Ex. smaragdina populations, sampled along a wide latitudinal range of Atlantic Forest, where the distribution of both species is co-occur; (ii) to make inferences on the evolutionary history of both species across the Atlantic Forest, and (iii) lastly, considering the wide distribution of both species through the Brazilian territory, we also investigated the evolutionary history of El. nigrita and Ex. smaragdina in a wider scenario, which included samples from others Brazilian forest domains.

MATERIALS AND METHODS

Study area and samplings

Bee samplings were carried out between 2012 and 2019, in 19 localities within the Atlantic Forest (Supplementary Material - Figure S1). Of these 19 localities, males of Ex. smaragdina and El. nigrita were collected in 12, not necessarily coincident (Supplementary Material - Table SI). As mentioned above, aiming to compare the set of haplotypes surveyed in areas of Atlantic Forest with other regions and biomes were both species are found, we included localities in the Amazon Forest (AM), and the Brazilian Caatinga and Cerrado biomes. One area of the Andean Montane Forest located in the municipality of Santa Fé de Antioquia in Colombia was included (Table SI). Two methods of collection were used, active collection (e.g., Nemésio 2010) and passive collection (PET bottle traps) (e.g., Santos Júnior et al. 2014, Martins et al. 2018), both using aromatic compounds (1,8-cineole, eugenol, methyl benzoate, methyl trans-cinnamate, methyl salicylate, skatole, p-tolyl and vanillin) for attraction. The specimens collected were stored in alcohol in the freezer at -20°C, or directly in the freezer. Finally, to complete the sampling, 15 pinned (dry material) males of Ex. smaragdina, deposited in the collection of insects from the “Centro de Coleções Taxonômicas da Universidade Federal de Minas Gerais” – CCT-UFMG, were included in the analyses (Table SII).

The DNA was extracted from the thoracic musculature or the hind leg of the bees, using the phenol-chloroform method (Sambrook & Russel 2001). In total, 118 specimens of El. nigrita and 81 of Ex. smaragdina were surveyed in our study. However, of these 97 of El. nigrita and 74 of Ex. smaragdina, since not all these specimens produced high-quality sequences (see Table SIII and Table SIV). The DNA pellet was re-suspended in 40 µl of LOW TE buffer (10 mM Tris-HCl, pH 8.0, 1 mM EDTA) and 1 μl was used for spectrophotometer quantification (NanoDrop Thermo Scientific 2000). After this, the material was diluted in a solution at 20 ng/μl concentration.

DNA processing

Two mitochondrial markers, COI and 16S genes, were chosen for our analysis, since they have shown some variation in other studies of orchid bees, e.g., COI (Dick et al. 2004, Nemésio et al. 2013) and 16S (Frantine-Silva et al. 2017). In addition, we also analyzed a segment of the Opsin nuclear gene (Michel-Salzat et al. 2004).

Amplifications of all genes were performed in a final volume of 15 μl PCR reaction mix, including 0.1 μl Taq Polymerase Platinum – 5u/μl (Platinum, Invitrogen, Brazil), 0.6 μl MgCl₂ – 50mM, 0.6 μl primers – 25 μM, 1.5 μl reaction buffer – 10X, 1.2 μl dNTPs – 2.5mM, and 1 μl DNA – ~20 ng /μl. The process was conducted in the thermocycler using different annealing temperatures (48-57°C see below). The amplification cycles were: 1 cycle of denaturation at 95°C for 5 min; 37 cycles of denaturation at 95°C for 45 s, primer annealing temperature for 30 s and extension at 72°C for 90 s; and the last pass of final extension at 72°C for 10 minutes. For mitochondrial DNA, the COI (Hebert et al. 2004) and 16S genes (Cameron et al. 1992) were used. In the case of nuclear DNA, the Opsin gene was chosen (Mardulyn & Cameron 1999) (Table SV). All PCR products were visualized in 0.8% agarose gel.

The positive samples were purified through the polyethylene glycol 20% (PEG) method according to Santos Júnior et al. (2015). For the sequencing reactions, the same primers were applied as in the PCR reactions. After purification, the PCR products were sequenced through the ABI 3130x1 sequencer following the manufacturer’s guidelines (Applied Biosystems). All the fragments were sequenced in forward and reverse directions. Sequenced fragments of the nuclear and mitochondrial genes were checked in Seqscape software 2.6 v. (Applied Biosystems, Darmstadt).

Genetic data analysis

The alignment of DNA sequences was made with the aid of the online tool MAFFT v. 7.475 (Katoh et al. 2017). For the alignment of the 16S sequence, the option Q-INS-i (the secondary structure of RNA is considered) was selected. After this, the Gblocks program (Dereeper et al. 2008) was used to remove regions of ambiguous alignment of this mitochondrial gene. Other fragments were performed in the MAFFT default.

Genetic diversity and population structure

For both species, we estimated the genetic diversity parameters, including the number of haplotypes (h), haplotype diversity (Hd), and nucleotide diversity (π), using the set of samples from each locality with DnaSP software (Librado & Rozas 2009). Firstly, we estimated these parameters based on sequences of the three gene-segments amplified for the set of samples from all localities. In this latter analysis genetic diversity measures were obtained for different clusters of each species, which were defined by spatial delimitation analyses performed in R software, using the package Geneland, version 4.0.8 run (Guillot et al. 2005). These same genetic diversity parameters were calculated separately for El. nigrita and Ex. smaragdina surveyed only in the Atlantic Forest areas.

Due to sample limitations, the uncorrelated haplotype frequency model (UHF) was chosen for some localities. The UHF model uses the Metropolis-Hasting algorithm to start from arbitrary values for all unknown parameters and to modify them in such a way that after many iterations, these values are close to the true values (Guillot et al. 2005). To identify genetic spatial discontinuities among populations of both species, the Geneland was used in three independent rounds. In the first, the number of possible population ranges from 1 to 10 was defined in the software. After this, based on the results, the software (i.e., K = 2) was executed again two times, using the number of clusters obtained for the program, to confirm the clustering. In all rounds, 10 million iterations were used of the Markov chain Monte Carlo (MCMC) and thinning by 100 to estimate the number of populations and geographic limits of the individuals.

The Analysis of Molecular Variance (AMOVA), calculated through 10,000 iterations, was carried out based only in sequences of COI gene, since amplified of this gene showed the highest genetic diversity for both orchid bee species. Based on the result from Geneland analysis, a hierarchical AMOVA was run, using the Tamura and Nei model, to consider variable base frequencies with equal transversion rates and variable transition rates, calculated using Arlequin for the results of each species (Excoffier & Lischer 2010).

To view the spatial distribution of the haplotypes that constructed the network, DnaSP was used to create the haplotype data file, removing the invariable sites to run in the Network v. 5 based on the Median-joining algorithm (Bandelt et al. 1999). This network was used to view the spatial distribution of haplotypes, e.g., the clusters found by Geneland.

The population structure from the species was examined based on the Φ_ST statistic, calculated in Arlequin. To test the existence of the correlation between genetic and geographical distances of observations across a landscape, the Mantel test in Alleles in Space software was used, with 10,000 permutations (Miller 2005).

Population demography and migration

The neutrality tests of Tajima’s D and Fu and Li’s D were carried out in Arlequin v. 3.5.2.2 (Excoffier & Lischer 2010). In both cases, only the COI sequences used were used in the analyses of both species.

The PartitionFinder 2 program was used to find the best partitioning scheme for each data set (Lanfear et al. 2016). From the PartitionFinder 2 results, Beast v. 2.6.0 software (Bouckaert et al. 2014) was run, using the Bayesian Skyline Plot method (BSP) to infer possible oscillations in the effective population size (Ne) over time between both species. Three independent runs were performed in Beast, one for each type of posterior distribution of data (normal, exponential, and strict). Each run was performed according to the following parameters for the dataset of each species: 100 million steps of MCMC, incrementing 100 steps, UPGMA initial tree, and F81+I (El. nigrita) and F81 (Ex. smaragdina) substitution model. To calibrate the molecular clock, the date of the most recent common ancestor of El. nigrita and Ex. smaragdina (9 mya; sd = 1 mya) was used, as shown only by Ramírez et al. (2010).

Tracer software was used to plot the BSP results of the posterior distribution (only when the effective sample size was ESS > 200) (Rambaut et al. 2013). The determination of best models of the posterior distribution of both species was based on the corrected Akaike Information Criterion (AICc).

To estimate the migration rate (M) between the clusters found in the Geneland results for each species, migration analysis was run in Migrate software, version 4.4.3 (Beerli 2016). This software uses the Metropolis-Hastings algorithm to calculate effective population sizes (Θ) and migration rates based on the coalescence theory. The following parameters were used for MCMC: inheritance scalers of 0.25, recording 50,000 steps, incrementing 1,000 steps, number of concurrent chains three replicates, visiting parameter values 150,000,000, and the first 10,000 genealogies were discarded as burn-in. The UPGMA tree was used as the starting genealogy and static heating with 4 automatic changes (temperatures were 1.0; 1.5; 3.0; 10¹¹) for each species.

Five models were evaluated to understand the migration dynamics of the groups found in Geneland (see Figures 1 and 3): (1) a model with one population (panmictic population); (2) – a full model with two population sizes and two migration rates, from cluster 1 to cluster 2, and from cluster 2 to cluster 1 (El. nigrita: AF+AM ↔ TN; Ex. smaragdina AF ↔ AM); (3) a model with two population sizes and one migration rate from cluster 2 to cluster 1 (El. nigrita: AF+AM ← TN; Ex. smaragdina (AF) ← (AM); (4) a model with two population sizes and one migration rate from cluster 1 to cluster 2 (El. nigrita: AF+AM → TN; Ex. smaragdina AF → AM; (5) a model with two separate groups of populations or independent evolutionary lineages. To find the best migration model indicated by Migrate for these species, the value of model probability from the Bezier approximation score was calculated for the different models (ΔBAS), as proposed by Beerli & Felsenstein (2001).

Figure 1
Population Clusters (K=2) identified using the UHF model implemented in Geneland Software for the COI gene of Eulaema nigrita. Cluster 1 includes the Amazon and Atlantic forests together, Cluster 2 showed areas of the Maranhão and Ceará in the area called, in this paper, Transitional area in northeastern Brazil (a), (b), and (c). The Lighter Colors Indicate Higher Probability Values >0.09. Diagrams and graph of the cluster formations (d).

Figure 3
Population Clusters (K=2) identified using the UHF model implemented in Geneland Software for the COI gene of Exaerete smaragdina. Cluster 1 includes the localities in the Amazon biome, while Cluster 2 showed localities in the Atlantic Forest (a), (b), and (c). The Lighter Colors Indicate Higher Probability Values >0.09. Diagrams and graph of the cluster formations (d).

RESULTS

Genetic diversity and population structure

In total, 97 males of El. nigrita were analyzed for COI (586 bp), 77 for 16S (~576 bp), and 58 for Opsin (585 bp) genes (Table I). The COI was conspicuously the most variable region for the number of haplotypes (h = 29) and haplotype diversity (Hd = 0.732). The highest value of nucleotide diversity was found for the 16S segment (π = 0.037). The neutrality test of the COI gene showed no significant values for El. nigrita (Table I).

Table I
Genetic diversity measures obtained for Eulaema nigrita and Exaerete smaragdina along the Atlantic Forest (above), and for samples of different 37 localities distributed in Brazilian territory and one site from Colombia (below). Diversity measures shown by samples of both species surveyed in the Atlantic Forest neutrality tests of Tajima’s D and Fu and Li’s D only for the most variable gene (COI). Sample sizes (N), number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), and standard deviation (sd). *p < 0.05; **p<0.01; ***p<0.001. Values for the clusters were generated in Geneland software. In the superior part of the table are shown the genetic diversity values obtained based only in COI sequences from samples of both species surveyed across Atlantic Forest. In the inferior part are shown results from mtDNA genes (16S, COI) and one nuDNA (Opsin) estimated from samples of 38 localities indicated in Table SI.

Species	COI gene	Genetic measures	AF only
Eulaema nigrita	COI (586 bp)	N	56
		h	12
		Hd (± sd)	0.438 (± 0.007)
		π	0.002
		Tajima’s D	-2.600***
		Fu and Li’s D	-4.642**
Exaerete smaragdina	COI (660 bp)	N	36
		h	8
		Hd (± sd)	0.556 (± 0.008)
		π	0.001
		Tajima’s D	-2.029*
		Fu and Li’s D	-3.089*
Species	Three genes	Genetic measures	Cluster 1	Cluster 2	Total
			AF + AM	TN
		N	64	33	97
Eulaema nigrita	COI (586 bp)	h	14	16	29
		Hd (± sd)	0.484 (± 0.006)	0.873 (± 0.002)	0.732 (± 0.002)
		π	0.007	0.012	0.017
		Tajima’s D	-1.411	0.179	0.840
		Fu and Li’s D	0.391	0.467	-0.763
	16S (~ 576 bp)	N	45	32	77
		h	3	2	3
		Hd (± sd)	0.279 (± 0.006)	0.125 (± 0.005)	0.234 (± 0.003)
		π	0.039	0.024	0.037
	Opsin (585 bp)	N	30	28	58
		h	8	5	9
		Hd (± sd)	0.595 (± 0.01)	0.381 (± 0.010)	0.495 (± 0.006)
		π	0.001	0.001	0.001
	Genes	Genetic measures	Cluster 1	Cluster 2	Total
			AF	AM
Exaerete smaragdina		N	36	38	74
	COI (660 bp)	h	9	7	14
		Hd (± sd)	0.551 (± 0.009)	0.485 (± 0.009)	0.677 (± 0.001)
		π	0.001	0.007	0.002
		Tajima’s D	-2.100*	-1.840*	-2.183*
		Fu and Li’s D	-3.311**	-1.541	-3.844**
	16S (~511 bp)	N	9	30	39
		h	2	6	8
		Hd (± sd)	0.222 (± 0.027)	0.310 (± 0.011)	0.567 (± 0.006)
		π	0.0004	0.001	0.001
	Opsin (668 bp)	N	3	16	19
		h	1	2	3
		Hd (± sd) (± sd)	0.000 (± 0.000)	0.400 (± 0.012)	0.433 (± 0.013)
		π	0.001	0.001	0.001

The spatial delimitation of the clusters based on the Geneland UHF model (only for COI gene) pointed to the existence of a range from K = 2 clusters for El. nigrita, with the haplogroup AF+AM encompassing all localities in the Amazon and Atlantic forests (AC2, AC4, AL, BA1, BA2, BA3, BA4, BA5, ES1, ES2, ES3, GO, MG1, MG2, MG3, MG4, and RO1) and only the Transitional areas (TN) in the northeast Brazil haplogroup in the states of Ceará and Maranhão (CE, MA1, MA2, MA3, MA4, MA5, and MA6) (Figure 1). The AF+AM cluster showed Hd = 0.484 and π = 0.0071. Analyzing separately the set of 56 males of El. nigrita sampled only along the AF, the number of haplotypes (h = 12) and estimates of Hd (0.438 ± 0.007) and π (0.002) were very close to those found for AF+AM. Distinctively from the results found for AF+AM, the neutrality values of Tajima’s (D = -2.029; p < 0.05) and Fu and Li’s tests (D = -3.089; p < 0.05), obtained for El. nigrita samples from AF, were statistically significant. Concerning to TN cluster, the value of Hd (0.873) found for TN was noticeably higher when compared to AF+AM cluster. Respectively, 13 and 10 private haplotypes were found to TN and AF+AM.

Among the total of 29 haplotypes found for the COI gene of El. nigrita, H1 (50.51%) was the most common among the localities sampled, being predominant in populations that inhabit the Amazon and Atlantic forests (Figure 2). On the other hand, 23 private haplotypes were identified for this species (Figure 2a). The Median-joining network haplotypes showed a starlike configuration close to H1, with several departing private haplotypes, evolving from one to two mutational steps, following a complex network with two medium vectors and several mutational steps (Figure 2b). The pattern in the frequency of H1 in these analyzed populations showed a breakup in the frequency of this haplotype between the humid tropical forests.

Figure 2
Twelve sites and frequency of haplotypes present in a fragment of the COI gene of Eulaema nigrita from Atlantic Forest and other different forest domains in the South America. (a) Pie charts indicate the frequency and distribution of each haplotype. Common haplotypes (black), haplotypes shared by at least 2 sites (gray), and private haplotypes (white). (b) Median-joining haplotype network for 29 haplotypes (H1-H29) for two haplogroups (Cluster 1 and 2), the blue color represents cluster 1 that includes the Amazon and Atlantic forests, while the yellow color is cluster 2 which represents the localities of northeast Brazil (CE, MA1, MA2, MA3, MA4, MA5, and MA6), and medium vectors (mv). Site codes are the same as in Table I, where AC, AL, BA, CE, ES, GO, MA, MG, and RO correspond to the usual abbreviations of the following Brazilian states Acre, Alagoas, Bahia, Ceará, Espírito Santo, Goiás. Maranhão, Minas Gerais, and Rondônia.

For Ex. smaragdina, 74 males were analyzed for the COI gene (660 bp), 39 for 16S (~511 bp), and 19 for Opsin (668 bp). In this species, the number of haplotypes and Hd were more variable for COI (h =14 haplotypes; Hd = 0.677), followed by 16S (h = 9; Hd = 0.567), and Opsin (h = 3; Hd = 0.433); the highest measure of nucleotide diversity (π = 0.002) was also found for the COI region (Table I). The demographic scenario analysis, performed through the neutrality tests of Tajima’s D and Fu and Li’s D, indicated population expansion only for Ex. smaragdina (COI: Tajima’s D= -2.1840, p < 0.01 and Fu and Li’s D = -3.8445, p < 0.01). In this species, the UHF identified two haplogroups (K = 2), which separate the Amazon and Atlantic Forest regions (Figure 3). The AF haplogroup includes the localities CE, PB1, PB2, PE1, PE2, PE3, PE4, BA1, BA2, ES1, and ES2 in the Atlantic Forest. The AM haplogroup includes localities in the Amazon Forest (AC1, AC2, AC4, AC5, CO, MA4, MA5, MS, PA, RO1, RO2, and SP). AM showed h = 7, Hd = 0.485, and π = 0.0079 while the AF haplogroup showed h= 9, Hd =0.551, and π = 0.0013. The private haplotypes for the two clusters were, respectively, seven and five haplotypes.

Analyzing the haplotype network for Ex. smaragdina, it is possible to notice that the H2 (41.89%) and H3 (39.18%) haplotypes were the most frequent (Figure 4a). The network haplotypes showed a double starlike configuration, with the H2 haplotype predominant in localities in the Amazon Forest and some sites in the Atlantic Forest (Figure 4b), while the H3 haplotype was present mainly in the population from the east and northeast of Brazil in the Atlantic Forest. This haplotype was also present in a single locality in west Brazil (RO2). Twelve private haplotypes from one to five mutational steps from H2 and H3 were randomly present along the sampling areas.

Figure 4
Twelve sample sites and frequency of haplotypes present in a fragment of the COI gene of Exaerete smaragdina from Atlantic Forest and other different forest domains in the South America. (a) Pie charts indicate the frequency and distribution of each haplotype. Haplotypes common (black), haplotypes shared by at least 2 sites (gray) and private haplotypes (white). (b) Median-joining network haplotype for 14 haplotypes (H1-H14) for two haplogroups (Cluster 1 and 2), light blue color represents the cluster 1 which is part the Atlantic forests with the localities CE, PB1, PB2, PE1, PE2, PE3, PE4, BA1, BA2, ES1, and ES2. Orange color in the network is represented by the other west of Brazil in the Amazon domain. Site codes are the same as in the Table I AC, BA, CE, ES, GO, MA, MS, PA, PB, PE, RO and SP correspond to usual abbreviations of the following Brazilian states, Acre, Bahia, Ceará, Espirito Santo, Maranhão, Mato Grosso do Sul, Pará, Paraíba, Pernambuco and Rondônia. The CO site correspond a localition of Santa Fe de Antioquia in the Colombia country.

For both species, the highest value of genetic structuring was found among populations within clusters (El. nigrita - Φ_ST = 0.840; p < 0.001 and Ex. smaragdina - Φ_SC = 0.374; p < 0.001) (Tables SVI and SVII). Our data point to different patterns of population structure in population levels for both species. In the case of the El. nigrita, the highest difference was found for the TN haplogroup (Table SVIII), whereas for the Ex. smaragdina the highest structuring was revealed for localities in the Amazon Forest compared to other points in the Atlantic Forest (Table SIX). Eulaema nigrita did not show a significant correlation between genetic and geographic distances (r = 0.031; p > 0.05). Exaerete smaragdina exhibited a positive correlation between these variables (Mantel test: r = 0.452; p < 0.05) (Figure S2), indicating the existence of greater gene flow among geographically closer populations.

Population demography and migration

The BSP revealed stability in effective population size (Ne) for both species in most of the Plio-Pleistocene period, with demographic expansion of about 0.5 mya (Figure 5). The full (AF+AM ↔ TN) model was selected by Migrate as the best model for El. nigrita, indicating that migration between the AF+AM and TN groups of populations was the more likely explanation for our data set. In the case of Ex. smaragdina, the best model shows that migration from the AM to AF (AF ← AM) groups of populations is more explanatory than the other hypotheses (Table II).

Figure 5
Bayesian Skyline Plot (BSP) based on the changes in effective population size of Eulaema nigrita (a) and Exaerete smaragdina (b). The strict clock was used to infer the demographic history of populations. The dark blue horizontal line shows median BSP estimate, and the blue area shows upper and lower 95% limits of the posterior density.

Table II
Comparison of migration models generated in Migrate 4.4.3 from the COI gene of two species of orchid bees (El. nigrita and Ex. smaragdina) in a cleptoparasite and host relationship. The arrows indicate the direction of migration between the groups, and * the best model.

Eulaema nigrita (AM+AF/ TN)
Models (n)	Bezier approximation score	ΔBAS	Model probability
1 Panmictic population	-1531,668.691	358,22627	0.000
2 (AF+AM) ↔ (TN)	-1352,555.556	0	0.994*
3 (AF+AM) ← (TN)	-1373,405.241	41,699.37	0.000
4 (AF+AM) → (TN)	-1357,654.873	10,198.634	0.006
Exaerete smaragdina (AM/ AF)
Models (n)	Bezier approximation score	ΔBAS	Model probability
1 Panmictic population	-1106,657.717	145,91671	0.000
2 (AF) ↔ (AM)	-1035,972.831	4,546.938	0.081
3 (AF) ← (AM)	-1033,699.362	0	0.786*
4 (AF) → (AM)	-1035,475.257	3,551.79	0.133

• *p < 0.05.

DISCUSSION

Pleistocene climate oscillations may have contributed to the current genetic structuring of the species studied here, as reported for other orchid bees (López-Uribe et al. 2014, Frantine-Silva et al. 2017), as well as other bee groups (Carvalho & Del Lama 2015, Miranda et al. 2016, 2017). The network configuration patterns, spatial delimitation reported by Geneland, migration models, and BSP, indicated that the distribution of genetic diversity in populations of El. nigrita and Ex. smaragdina throughout the Atlantic Forest would have a distinct evolutionary history relative to the TN and Amazon Forest, respectively. In both cases, possibly influenced by isolations and subsequent population expansions from refugia areas, enabling demographic expansion of both species during the Pleistocene, and corroborating with data reported for bees (Santos Júnior et al. 2015, Miranda et al. 2016, 2017), and other organisms (Haffer 1969, Carnaval et al. 2009).

Genetic structure and genetic diversity

Although the present study focused on information from the mitochondrial COI region, which is less variable than non-genic regions (introns and microsatellites), the haplotypic diversity values found for the TN cluster of El. nigrita were higher when compared to what has been reported for the mitochondrial DNA of other orchid bee species (López-Uribe et al. 2014, Frantine-Silva et al. 2017). At the same time, data obtained here for both species also suggest clues to recent population expansions throughout the Atlantic Forest, markedly influencing the distribution of genetic diversity among populations of El. nigrita and Ex. smaragdina.

Population demography and migration

Data from the present study place the population expansions of Ex. smaragdina and El. nigrita during the Pleistocene, a period that included sequences of glaciation events that changed forest patterns around the world. The cycles of glaciations promoted expansion and retraction of the arid environments of South America (Haffer 1969, Werneck 2011), providing opportunities for species from more forested environments (e.g., Ex. smaragdina) and adapted to open areas (e.g., El. nigrita) to establish (Carvalho & Almeida 2011).

Dick et al. (2004), studying the possible effect of the Andean cordilleras on the same species, did not find a phylogeographical structure in Amazon euglossines. This fact may be explained by the recent population expansion of the bees through different environments in South America. Another possible explanation for this expansion may be the permanence of large areas of Amazon formation that remained unchanged by the events of climatic oscillations in the Pleistocene (Colinvaux et al. 2000). Following this line, the BSP showed stability between 2 and 7 mya for the two species (see Figure 5), which may have represented a period of great climatic instability on the demographic pattern of populations during the Plio-Pleistocene. There was expansion in the effective size of both species only in the late Pleistocene, when the climatic changes seem to have similarly affected the historical demography patterns of both El. nigrita and Ex. smaragdina, as reported for other bee species (Miranda et al. 2016, 2017, Frantine-Silva et al. 2017). When there was climate stability, a higher Ne (Ne=~80) was observed for the Ex. smaragdina compared to El. nigrita. Concerning this fact, it has been suggested that El. nigrita is not the only host of Ex. smaragdina throughout the distribution of both species (Nemésio & Silveira 2006). These latter authors highlight for the high abundance of Ex. smaragdina in several areas where El. nigrita is absent in the Amazon Basin, pointing out El. cingulata (Fabricius) or El. meriana (Olivier) as potential hosts of Ex. smaragdina, a supposition still to be proved (Nemésio & Silveira 2006).

Overall, clearer evidence of population expansion was obtained for Ex. smaragdina, which revealed significant values for Tajima’s D and Fu and Li’s D, as well as the starlike configuration in the haplotype network. In addition, the combination Hd > 0.5 together with π < 0.5%, according to Grant & Bowen (1998), may suggest a population bottleneck followed by rapid population growth and accumulation of mutations. In the case of El. nigrita, the population expansion was also indicated by both Tajima’s D and Fu and Li’s D estimates and the starlike configuration, readily apparent among the Atlantic Forest in the haplotype network, including most haplotypes separated by only one mutational step. On the other hand, the haplotypic (Hd < 0.5) and nucleotide (π > 0.5%) diversities of the AF+AM cluster suggest divergence among geographically subdivided populations (Grant & Bowen 1998), which seems plausible since the AF+AM analysis included haplotypes H3 (RO1) and H2 (GO) that belong to the TN cluster haplogroup, and this possibly influenced the genetic diversity estimates and neutrality tests. In any case, BSP results agreed that both species possibly had expansion in the effective size in the late Pleistocene, starting about 500,000 years ago, although the times of onset of expansion and the routes taken may present important differences between the species.

In the case of El. nigrita, comparisons between AF+AM and TN brought information within the evolutionary history of the species. The TN group, besides showing haplotypic diversity twice as high as AF+AM, and with 13 of its 16 haplotypes being private, also showed no signs of recent evolutionarily population expansions. In fact, the Hd > 0.5 and π > 0.5% values obtained for TN suggest a large and stable population with a long evolutionary history or secondary contact between different lineages (Grant & Bowen 1998). At the same time, TN showed no significant values in neutrality tests and its haplotypes are less related than those in AF+AM, including several mutational steps, indicating distinct and longer-established populations than those in the Atlantic Forest. Due to failures during the sequencing process of the Ex. smaragdina samples, the sampling of the TN region for this species was very limited, which affects comparisons for this area.

Indeed, expansions from populations that remained in refuge areas during the Pleistocene climatic oscillations seem plausible for both species in the present study. This scenario is not new and has already been reported for stingless bees (Carvalho & Del Lama 2015, Miranda et al. 2016, 2017), as well as in other taxa, such as mammals (Vivo 1997, Leite et al. 2016), lizards (Werneck 2011, Pellegrino et al. 2011), birds (Haffer 1969, Batalha-Filho et al. 2013), and plants (Pennington et al. 2004, Martini et al. 2007), reinforcing the role of these areas as a shelter for high genetic diversity.

The Pleistocene models of identification of stable (refuge areas) versus unstable areas are one of the hypotheses to explain the high diversity of species in the Neotropical region (Haffer 1969, Carnaval & Moritz 2008, Carnaval et al. 2009). These areas have an important role in the conservation of bees, sheltering high genetic diversity (Carvalho & Del Lama 2015). In this sense, some studies on orchid bees in different scales such as transcontinental (López-Uribe et al. 2014) and Atlantic Forest areas found these areas to be important refuges for orchid bee species (Garraffoni et al. 2017, Miranda et al. 2019). Among the other several theories trying to explain the biodiversity of South America, the Forest refuge hypothesis is the most accepted (Haffer 1969).

In the Pleistocene, refuges corresponded to stable environments during the climatic change of this epoch (Carnaval & Moritz 2008, Carnaval et al. 2009). López-Uribe et al. (2014) tested the effect of climatic instability of the Pleistocene in three Eulaema species, using mitochondrial (COI and Cyt b) and nuclear (microsatellites) markers allied to niche modeling. The authors suggested three predictions for the orchid bee species with respect to historical climatic instability; (1) there was an uneven reduction in the number of forest refuge areas during the dry period of the Pleistocene for all orchid bee species, (2) orchid bees species with low physiological tolerance may show a stronger spatial phylogeographical structure following the refuge hypothesis, and (3) the spatial patterns of genetic diversity of these bees are the consequence of isolation and colonization events during several cycles of forest contraction and expansion. Both species studied here seem to fit these hypotheses.

Although the genetic data from the present study do not clearly show the refuge area from which the population of El. nigrita expanded in the Atlantic Forest, it seems quite plausible that TN is within or near one of these areas, maintaining ancestral genetic diversity patterns and possibly influencing other areas in the Amazon region, as suggested by haplotype H3 in RO1 and H2 in GO. In fact, López-Uribe et al. (2014) found a suitable area during the Last Glacial Maxima (LGM, 21kya) on the Maranhão coast (TN cluster) for El. cingulata. Currently, El. nigrita and El. cingulata cohabit in several ecosystems from rainforest environments to open formations of Cerrado and Caatinga in South America (Silveira et al. 2015, Martins et al. 2018), which may suggest possible cohabitation of these species in the suitable area of the Maranhão coast during the LGM, in the same area as the TN cluster (López-Uribe et al. 2014). Studying refuge areas for stingless bees (Meliponini), Camargo & Pedro (2003) also pointed out areas in Maranhão and Pará states as regions of climate stability during the Pleistocene. In addition, Martins et al. (2021) found high levels of genetic structuring for Euglossa cordata (Linnaeus, 1758), reinforcing the importance of the TN refuge in this region for bee conservation.

Starting from the idea that the Amazon Basin is the original center of the Euglossini (Dressler 1982, Ramírez et al. 2010), it is possible that El. nigrita was isolated in the TN due to unfavorable environmental conditions during the glaciation events in the Pleistocene. After these climate oscillations, the species studied found conditions to establish themselves and exhibited high demographic expansion due to the high dispersion capacity of the males, in different moments during the Pleistocene. Concurrently, it seems clear that populations of El. nigrita from the Atlantic Forest expanded from a refuge area different from TN, including an evolutionary history that involves the previous dispersion of the species to the east coast of South America, as demonstrated by Miranda et al. (2019) who identified priority areas for the conservation of ten orchid bee species, among these the two studied here, through ecological niche modeling.

In the case of Ex. smaragdina, genetic data suggest expansion from at least two refuges, one influencing the Amazon populations more markedly and the other more evident among the Atlantic Forest populations. Although more conclusive discussions require further studies, one of the areas of highest genetic diversity for Ex. smaragdina was in one of the Pernambuco sites (PE3), precisely within one of the suggested refuge areas for the Atlantic Forest suggested for both species by Miranda et al. (2019). In this context, it seems plausible that the current distribution of genetic diversity of mitochondrial DNA of Ex. smaragdina included reflections of an expansion from Pernambuco, dispersing over generations towards the south of the Atlantic Forest formation. Concurrently, populations still dispersing in the north-south direction could be facing more difficulty establishing in the current scenario, and this may be contributing to the lower frequency in smaller latitudes.

In some Euglossini genera, such as Eulaema and Exaerete, the dispersion capacity is directly related to body size (Janzen 1971, Dick et al. 2004, Pokorny et al. 2015). In this case, larger bees have higher advantages such as thermal tolerance, high energy to fly long-distances, and the capacity to colonize new areas (founder effect) (Dick et al. 2004). Thus, the best migration model was full (AF+AM ↔ TN) for El. nigrita, indicating that migration between clusters occurs in both directions. For Ex. smaragdina, the best model showed dispersion in the direction (AF ← AM), and this pattern appears to reflect the historical pattern that the whole group follows (Nemésio 2009), which seems plausible since the Amazon is reported as the original center of the Euglossini (Dressler 1982, Ramírez et al. 2010). Eulaema nigrita did not show this pattern because AM and AF were grouped, and TN showed different haplotypes.

In fact, AM and AF clusters examined for Ex. smaragdina showed similar values of genetic diversity and very close haplogroups separated by only one mutation. This also indicates common founder populations for AM and AF, and these data suggest that the populations of this species are well established in these forest formations (Nemésio & Silveira 2006). However, the lack of higher genetic diversity values is probably the consequence of the non-sampling of the original center of Ex. smaragdina, due to limitations in our sampling and failures during the sample sequencing process.

According to López-Uribe et al. (2014), the higher structuration observed in mtDNA and correlation of geographical and genetic distance is the result of a different dispersion capacity of Euglossini males (up to 40 km in a few days) (Dressler 1982, Pokorny et al. 2015) and females (up to 23 km in the same day) (Janzen 1971). Moreover, the genetic transmission of the material of the mtDNA that is exclusively maternal (Avise 2009), associated with the higher intraspecific mutation rate, is the possible cause of the structuration reported for some studies (López-Uribe et al. 2014, Penha et al. 2015).

Interestingly, the samples of Ex. smaragdina in the present study showed a direct relationship between genetic and geographical variables (Mantel test: r = 0.452; p<0.05), corroborating isolation by distance, which is possibly influenced by the dispersal pattern after population expansion. A clear pattern of genetic divergence was observed for populations separated by more than 418 km, highlighting samples from Espírito Santo with the highest isolation levels. In the case of El. nigrita, the absence of this correlation was already expected, since TN and AF+AM populations include different haplogroups.

CONCLUSION

Despite the closely related ecological relationship between El. nigrita and Ex. smaragdina, our data reveal different patterns of structure for these orchid bee species, suggesting different evolutionary history for both species throughout their distribution. On the other hand, the concomitant patterns of population expansion for the studied species, such as the starlike configuration, the presence of few mutations between haplotypes, and significant values for neutrality tests suggest similar response patterns to climate oscillations during the Pleistocene for both species. The highest values of genetic diversity detected for El. nigrita within the Atlantic Forest (ES1 and ESP3) and TN region, outside the AF, point out these areas as potential refuges for this during Pleistocene. For Ex. smaragdina, Pernambuco (PE3) may have provided suitable conditions for the survival and maintenance of the high genetic diversity of orchid bees, as highlighted for other orchid bee species. Thus, the preservation of these environments becomes vital for the maintenance of the populations of these important elements of the Neotropical fauna.

ACKNOWLEDGMENTS

We are grateful for the financial support of the Fundação de Amparo à Pesquisa e ao Desenvolvimento Científico e Tecnológico do Maranhão (FAPEMA Edital Nº 031/ 2016) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES - Finance Code 001). DCM thanks to CAPES the scholarship awarded; SHS receives a productivity research fellowship from Conselho Nacional de Desenvolvimento Científico e Tecnológico - CNPq (PQ 305343/2018-1). We thank the collaborator at the Laboratory of Animal Genetics and Ecology (LAGEA/UEL) Thais Kotelok-Diniz for the support in genetic analysis. We also thank the anonymous reviewers for their valuable comments.

SUPPLEMENTARY MATERIAL

Tables SI, SII, SIII, SIV, SV, SVI, SVII, SVIII, SIX.

REFERENCES

• AGUIAR WM & GAGLIANONE MC. 2012. Euglossine bee communities in small forest fragments of the Atlantic Forest, Rio de Janeiro state, southeastern Brazil (Hymenoptera, Apidae). Rev Bras Entomol 56(2): 210-219.
• AVISE JC. 2009. Phylogeography: retrospect and prospect. J Biogeogr 36: 3-15.
• BANDELT HJ, FORSTER P & ROHL A. 1999. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol 16: 37-48.
• BATALHA-FILHO H, FJELDSÅ J, FABRE PH & MIYAKI CY. 2013. Connections between the Atlantic and the Amazonian Forest avifaunas represent distinct historical events. J Ornithol 154: 41-50.
• BATALHA-FILHO H, WALDSCHMIDT AM, CAMPOS LAO, TAVARES MG & FERNANDES-SALOMÃO TM. 2010. Phylogeography and historical demography of the Neotropical stingless bee Melipona quadrifasciata (Hymenoptera, Apidae): incongruence between morphology and mitochondrial DNA. Apidologie 41: 534-547.
• BEERLI P. 2016. MIGRATE-N. v. 4.4.3. Computer program and documentation distributed by the author. Available in: https://peterbeerli.com/migrate-html5/index.html
  » https://peterbeerli.com/migrate-html5/index.html
• BEERLI P & FELSENSTEIN J. 2001. Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach. Proc Natl Acad Sci USA 98: 4563-4568.
• BOUCKAERT R, HELED J, KÜHNERT D, VAUGHAN T, WU C-H, XIE D, SUCHARD MA, RAMBAUT A & DRUMMOND AJ. 2014. BEAST 2: A Software Platform for Bayesian Evolutionary Analysis. PLoS Comput Biol 10: e1003537.
• CAMARGO JMF & PEDRO SRM. 2003. Meliponini neotropicais: o gênero Partamona Schwarz, 1939 (Hymenoptera, Apidae, Apinae) - bionomia e biogeografia. Rev Bras Entomol 47: 311-372.
• CAMERON SA, AUSTIN AD & WHARTON RA. 1992. The application of nucleotide sequence data to phylogeny of the Hymenoptera: a review. J Hymenopt Res 1: 63-79.
• CARNAVAL AC, HICKERSON MJ, HADDAD CFB, RODRIGUES MT & MORITZ C. 2009. Stability predicts genetic diversity in the Brazilian Atlantic Forest hotspot. Science 323: 785-789.
• CARNAVAL AC & MORITZ C. 2008. Historical climate modelling predicts patterns of current biodiversity in the Brazilian Atlantic Forest. J Biogeogr 35: 1187-1201.
• CARVALHO CJB & ALMEIDA EAB. 2011. Biogeografia da América do Sul: Padrões & Processos, 1st ed. São Paulo: Ed. Roca, 306 p.
• CARVALHO AF & DEL LAMA MA. 2015. Predicting priority areas for conservation from historical climate modelling: stingless bees from Atlantic Forest hotspot as a case study. J Insect Conserv 19: 581-587.
• COLINVAUX PA, DE OLIVEIRA PE & BUSH MB. 2000. Amazonian and neotropical plant communities on glacial time-scales: The failure of the aridity and refuge hypotheses. Quaternary Sci Rev 19: 141-169.
• CORDEIRO GD, BOFF S, CAETANO TA, FERNANDES PC & ALVES-DOS-SANTOS I. 2013. Euglossine bees (Apidae) in Atlantic forest areas of São Paulo State, southeastern Brazil. Apidologie 44: 254-267.
• DANFORTH BN, MICKLEY RL & NEFF JL. 2019. The Solitary Bees. Biology, Evolution, Conservation. Princeton: Princeton University Press, 472 p.
• DEREEPER A ET AL. 2008. Phylogeny.fr: robust phylogenetic analysis for the non-specialist. Nucl Acid Res 36: W465-W469.
• DICK CW, ROUBIK DW, GRUBER KF & BERMINGHAM E. 2004. Long-distance gene flow and cross-Andean dispersal of lowland rainforest bees (Apidae: Euglossini) revealed by comparative mitochondrial DNA phylogeography: phylogeography of euglossine bees. Mol Ecol 13: 3775-3785.
• DRESSLER RL. 1982. Biology of the orchid bees (Euglossini). Ann Rev Ecol Syst 13: 373-394.
• EXCOFFIER L & LISCHER HEL. 2010. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour 10: 564-567.
• FRANTINE-SILVA W, GIANGARELLI DC, PENHA RES, SUZUKI KM, DEC E, GAGLIANONE MC, ALVES-DOS-SANTOS I & SOFIA SH. 2017. Phylogeography and historical demography of the orchid bee Euglossa iopoecila: signs of vicariant events associated to Quaternary climatic changes. Conserv Genet 18: 539-552.
• GARÓFALO CA & ROZEN JG. 2001. Parasitic behavior of Exaerete smaragdina with descriptions of its mature oocyte and larval instars (Hymenoptera: Apidae: Euglossini). Am Mus Novit 3349: 1-28.
• GARRAFFONI ARS, MOURA FR & LOURENÇO AP. 2017. Areas of endemism in the Atlantic Forest: quantitative biogeography insights from orchid bees (Apidae: Euglossini). Apidologie 48: 513-522.
• GRANT WS & BOWEN BW. 1998. Shallow population histories in deep evolutionary lineages of marine fishes: insights from sardines and anchovies and lessons for conservation. J Hered 89: 415-426.
• GUILLOT G, MORTIER F & ESTOUP A. 2005. Geneland: a computer package for landscape genetics. Mol Ecol Notes 5: 712-715.
• HAFFER J. 1969. Speciation in Amazonian Forest birds. Science 165: 131-137.
• HEBERT PDN, PENTON EH, BURNS JM, JANZEN DH & HALLWACHS W. 2004. Ten species in one: DNA barcoding reveals cryptic species in the neotropical skipper butterfly Astraptes fulgerator. Proc Natl Acad Sci USA 101: 14812-14817.
• JANZEN DH. 1971. Euglossine Bees as Long-Distance Pollinators of Tropical Plants. Science 171: 203-205.
• KATOH K, ROZEWICKI J & YAMADA KD. 2017. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform 20: 1160-1166.
• LANFEAR R, FRANDSEN PB, WRIGHT AM, SENFELD T & BRETT C. 2016. PartitionFinder 2: New Methods for Selecting Partitioned Models of Evolution for Molecular and Morphological Phylogenetic Analyses. Mol Biol Evol 34: 772-773.
• LEITE YLR ET AL. 2016. Neotropical forest expansion during the last glacial period challenges refuge hypothesis. Proc Natl Acad Sci USA 113: 1008-1013.
• LIBRADO P & ROZAS J. 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451-1452.
• LOPES EJM ET AL. 2022. Local abundance of neotropical orchid bees in Amazon forests not related to large-scale climate suitability. Insect Conserv Divers 15: 693-703.
• LÓPEZ-URIBE MM, ZAMUDIO KR, CARDOSO CF & DANFORTH BN. 2014. Climate, physiological tolerance and sex-biased dispersal shape genetic structure of Neotropical orchid bees. Mol Ecol 23: 1874-1890.
• MARDULYN P & CAMERON SA. 1999. The Major Opsin in Bees (Insecta: Hymenoptera): A Promising Nuclear Gene for Higher Level Phylogenetics. Mol Phylogenet Evol 12: 168-176.
• MARTINI AMZ, FIASCHI P, AMORIM AM & PAIXÃO JL. 2007. A hot-point within a hot-spot: a high diversity site in Brazil’s Atlantic Forest. Biodivers Conserv 16: 3111-3128.
• MARTINS DC, ALBUQUERQUE PMC, SILVA FS & REBÊLO JMM. 2018. Orchid bees (Apidae: Euglossini) in Cerrado remnants in northeast Brazil. J Nat Hist 52: 627-644.
• MARTINS DC, ALBUQUERQUE PMC, REBÊLO JMM, KOTELOK-DINIZ T, SOFIA SH & FRANTINE-SILVA W. 2021. Phytogeographic regions and geographic distance do not predict genetic structure in the orchid bee Euglossa cordata. J Apic Res 1-12. https://doi.org/10.1080/00218839.2021.1905373
  » https://doi.org/10.1080/00218839.2021.1905373
• MICHEL-SALZAT A, CAMERON SA & OLIVEIRA ML. 2004. Phylogeny of the orchid bees (Hymenoptera: Apinae: Euglossini): DNA and morphology yield equivalent patterns. Mol Phylogenet Evol 32: 309-323.
• MILLER MP. 2005. Alleles in space (AIS): computer software for the joint analysis of interindividual spatial and genetic information. J Hered 96: 722-724.
• MIRANDA EA, BATALHA-FILHO H, CONGRAINS C, CARVALHO AF, FERREIRA KM & DEL LAMA MA. 2016. Phylogeography of Partamona rustica (Hymenoptera, Apidae), an endemic stingless bee from the Neotropical Dry Forest Diagonal. PLoS ONE 11: e0164441. https://doi.org/10.1371/journal.pone.0164441.
  » https://doi.org/10.1371/journal.pone.0164441
• MIRANDA EA, CARVALHO AF, JESUS GOMES-MIRANDA J, SOUZA CR & COSTA MA. 2019. Priority areas for conservation of orchid bees (Apidae, Euglossini) in the Atlantic Forest. J Insect Conserv 23: 613-621.
• MIRANDA EA, FERREIRA KM, CARVALHO AT, MARTINS CF, FERNANDES CR & DEL LAMA MA. 2017. Pleistocene climate changes shaped the population structure of Partamona seridoensis (Apidae, Meliponini), an endemic stingless bee from the Neotropical dry forest. PLoS ONE 12: e0175725.
• NEMÉSIO A. 2009. Orchid bees (Hymenoptera:Apidae) of the Brazilian Atlantic Forest. Zootaxa 2041: 1-242.
• NEMÉSIO A. 2010. The orchid-bee fauna (Hymenoptera: Apidae) of a forest remnant in northeastern Brazil, with new geographic records and an identification key to the known species of the Atlantic Forest of northeastern Brazil. Zootaxa 2656: 55-66.
• NEMÉSIO A & SANTOS JÚNIOR JE. 2013. Eufriesea zhangi sp. n. (Hymenoptera: Apidae: Euglossina), a new orchid bee from Brazil revealed by molecular and morphological characters. Zootaxa 3609: 568-582.
• NEMÉSIO A & SILVEIRA FA. 2006. Deriving ecological relationships from geographical correlations between host and parasitic species: an example with orchid bees. J Biogeogr 33: 91-97.
• PELLEGRINO KCM, RODRIGUES MT, HARRIS DJ, YONENAGA-YASSUDA Y & SITES JR JW. 2011. Molecular phylogeny, biogeography and insights into the origin of parthenogenesis in the Neotropical genus Leposoma (Squamata: Gymnophthalmidae): Ancient links between the Atlantic Forest and Amazonia. Mol Phylogenet Evol 61: 446-459.
• PENHA RES, GAGLIANONE MC, ALMEIDA FS, BOFF SV & SOFIA SH. 2015. Mitochondrial DNA of Euglossa iopoecila (Apidae, Euglossini) reveals two distinct lineages for this orchid bee species endemic to the Atlantic Forest. Apidologie 46: 346-358.
• PENNINGTON RT, LAVIN M, PRADO DE, PENDRY CA, PELL SK & BUTTERWORTH CA. 2004. Historical climate change and speciation: neotropical seasonally dry forest plants show patterns of both Tertiary and Quaternary diversification. Philos Trans R Soc Lond B Biol Sci 359: 515-538.
• POKORNY T, LOOSE D, DYKER G, QUEZADA-EUÁN JJG & ELTZ T. 2015. Dispersal ability of male orchid bees and direct evidence for long-range flights. Apidologie 46: 224-237.
• RAMBAUT A, SUCHARD MA & DRUMMOND AJ. 2013. Tracer v1.6: Available in: http://tree.bio.ed.ac.uk/software/TRACER/
  » http://tree.bio.ed.ac.uk/software/TRACER/
• RAMÍREZ SR, ROUBIK DW, SKOV C & PIERCE NE. 2010. Phylogeny, diversification patterns and historical biogeography of euglossine orchid bees (Hymenoptera: Apidae). Biol J Linn Soc 100: 552-572.
• REBÊLO JMM & GARÓFALO CA. 1997. Comunidades de machos de Euglossini (Hymenoptera: Apidae) em Matas Semidecíduas do Nordeste do Estado de São Paulo. An Soc Entomol Brasil 26(2): 243-255.
• ROUBIK DW & HANSON PE. 2004. Orchid bees from Tropical America. Biology and field guide. Santo Domingo de Heredia: INBio Press, 370 p.
• SAMBROOK J & RUSSEL DW. 2001. Molecular Cloning - A Laboratory Manual-Cold Spring Harbor Laboratory Press. Cold Spring Harbor, New York: CSH Laboratory Press, 2100 p.
• SANTOS JÚNIOR J, FERRARI R & NEMÉSIO A. 2014. The orchid-bee fauna (Hymenoptera: Apidae) of a forest remnant in the southern portion of the Brazilian Amazon. Braz J Biol 74: S184-S190.
• SANTOS JÚNIOR JE, SANTOS FR & SILVEIRA FA. 2015. Hitting an Unintended Target: Phylogeography of Bombus brasiliensis Lepeletier, 1836 and the first new Brazilian bumblebee species in a century (Hymenoptera: Apidae). PLoS ONE 10: e0125847.
• SILVEIRA GC, FREITAS RF, TOSTA THA, RABELO LS, GAGLIANONE MC & AUGUSTO SC. 2015. The orchid bee fauna in the Brazilian savanna: do forest formations contribute to higher species diversity? Apidologie 46: 197-208.
• STORCK-TONON D, MORATO EF, MELO AWF & OLIVEIRA ML. 2013. Orchid bees of forest fragments in Southwestern Amazonia. Biota Neotropica 13: 133-141.
• STORCK-TONON D, MORATO EF & OLIVEIRA ML. 2009. Fauna de Euglossina (Hymenoptera: Apidae) da Amazônia Sul-Ocidental, Acre, Brasil. Acta Amazon 39: 693-706.
• VIVO M. 1997. Mammalian evidence of historical ecological change in the Caatinga semiarid vegetation of northeastern Brazil. J Comput Biol 2: 65-73.
• WCISLO WT. 1987. The roles of seasonality, host synchrony, and behaviour in the evolution and distribution of nest parasites in Hymenoptera (Insecta), with special reference to bees (Apoidea). Biol Rev 62: 515-542.
• WERNECK FP. 2011. The diversification of eastern South American open vegetation biomes: Historical biogeography and perspectives. Quaternary Sci Rev 30: 1630-1648.
• ZUCCHI R, SAKAGAMI SF & CAMARGO JMF. 1969. Biological observations on a Neotropical parasocial bee, Eulaema nigrita, with a review on the biology of Euglossinae (Hymenoptera, Apidae). J New York Entomol Soc 17: 271-380.

Publication Dates

• Publication in this collection
  20 Oct 2023
• Date of issue
  2023

History

• Received
  8 July 2022
• Accepted
  23 Sept 2022