The orchid bees constitute a clade of prominent insect pollinators distributed throughout the Neotropical region. Males of all species collect fragrances from natural sources, including flowers, decaying vegetation and fungi, and store them in specialized leg pockets to later expose during courtship display. In addition, orchid bees provide pollination services to a diverse array of Neotropical angiosperms when foraging for food and nesting materials. However, despite their ecological importance, little is known about the evolutionary history of orchid bees. Here, we present a comprehensive molecular phylogenetic analysis based on ∼4.0 kb of DNA from four loci [cytochrome oxidase (CO1), elongation factor 1-α (EF1-α), arginine kinase (ArgK) and RNA polymerase II (Pol-II)] across the entire tribe Euglossini, including all five genera, eight subgenera and 126 of the approximately 200 known species. We investigated lineage diversification using fossil-calibrated molecular clocks and the evolution of morphological traits using disparity-through-time plots. In addition, we inferred past biogeographical events by implementing model-based likelihood methods. Our dataset supports a new view on generic relationships and indicates that the cleptoparasitic genus Exaerete is sister to the remaining orchid bee genera. Our divergence time estimates indicate that extant orchid bee lineages shared a most recent common ancestor at 27–42 Mya. In addition, our analysis of morphology shows that tongue length and body size experienced rapid disparity bursts that coincide with the origin of diverse genera (Euglossa and Eufriesea). Finally, our analysis of historical biogeography indicates that early diversification episodes shared a history on both sides of Mesoamerica, where orchid bees dispersed across the Caribbean, and through a Panamanian connection, thus reinforcing the hypothesis that recent geological events (e.g. the formation of the isthmus of Panama) contributed to the diversification of the rich Neotropical biota.
For more than 150 years, orchid bees (tribe Euglossini) have caught the attention of scientists and naturalists alike, mainly as a result of the intricate association of male bees as pollinators of orchid flowers (Dressler, 1967, 1968, 1982a; Roubik & Hanson, 2004). Male orchid bees collect aromatic compounds from orchids (and other nonfloral sources) and store them in specialized pouches located in their hind tibiae (Fig. 1). They later expose such fragrances by presenting small samples of scent to females during courtship (Eltz et al., 1999; Eltz, Roubik & Whitten, 2003; Bembé, 2004, Eltz, Sager & Lunau 2005a, Eltz, Roubik & Lunau, 2005b). Fragrances are thus thought to act as species-specific pheromone analogues that evolved in response to reproductive interference among closely related lineages (Zimmermann, Ramírez & Eltz, 2009). Whilst collecting fragrances, male orchid bees visit and pollinate a large diversity of orchids (Williams, 1982; Williams & Whitten, 1983), as well as other plant families (Ramírez, Dressler & Ospina, 2002). Since the 1960s, when fragrance-gathering behaviour was discovered (Vogel, 1963, 1966; Dodson et al., 1969), synthetic chemical baits that attract male bees have been widely used in ecological surveys across the Neotropics (Roubik & Ackerman, 1987; Roubik & Hanson, 2004). These surveys revealed a great deal about orchid bee ecology, population biology and natural history (Dressler, 1982a; Ramírez et al., 2002; Cameron, 2004; Roubik & Hanson, 2004; Ramírez, 2009).
The tribe Euglossini belongs to the bee family Apidae, and is contained within a clade commonly known as the corbiculate bees (Michener, 2007). Female corbiculate bees have basket-shaped hind tibiae (corbiculae) for transporting pollen and nesting materials to the nest (Roubik, 1989), a trait that has been secondarily lost in nest-parasites (Michener, 2007). In addition to orchid bees, corbiculate bees also include the primitively eusocial bumble bees (tribe Bombini) and the highly eusocial honey bees (tribe Apini) and stingless bees (tribe Meliponini) (Michener, 2007). However, unlike their close relatives, orchid bees do not display eusocial traits, such as worker sterility or cooperative brood care, and instead exhibit solitary, communal or parasocial nesting (Zucchi, Sakagami & Camargo, 1969; Cameron & Ramírez, 2001; Soucy, Giray & Roubik, 2003; Augusto & Garófalo, 2004; Roubik & Hanson, 2004).
Corbiculate bees are contained within the monophyletic long-tongued bees (Michener, 1990; Danforth et al., 2006), a status supported by both morphological and molecular characters (Michener, 2007; Kawakita et al., 2008), but the phylogenetic position of Euglossini within corbiculate bees has been a subject of considerable debate. Because all tribes of corbiculate bees, except Euglossini, exhibit some degree of advanced eusociality (Michener, 2007), the position of orchid bees with respect to the remaining corbiculate bee lineages can potentially establish the number of origins of eusocial behaviour in this group (for a most up-to-date review, see Kawakita et al. 2008). Whereas hypotheses based on morphological characters tend to favour orchid bees as sister to the rest of the corbiculate bees, and thus suggest a single origin of eusociality among corbiculate clades (Michener, 1944; Engel, 2001; but see Winston & Michener, 1977), most studies based on molecular characters favour the hypothesis of orchid bees as sister to honey bees, thus implying a dual origin of highly eusocial behaviour (Cameron & Mardulyn, 2001; Kawakita et al., 2008).
The Euglossini contains five (Fig. 1) well-delineated extant genera (Aglae, Eufriesea, Euglossa, Exaerete and Eulaema) and over 200 described species (Moure, Melo & Faria, 2007). Several alternative hypotheses regarding the relationships between orchid bee genera have been advanced on the basis of both molecular and morphological characters (Fig. 2A). Moreover, classification schemes and phylogenetic hypotheses from morphology have been proposed for portions of all genera.
The genus Aglae contains a single large species (A. caerulea), which is an exclusive cleptoparasite (or ‘cuckoo’ parasite) of Eulaema (Dressler, 1982a). Aglae caerulea is restricted to the lowland wet forests of the Chocó region and Amazon basin of South America (Ramírez et al., 2002; Anjos-Silva, Camilo & Garófalo, 2006). The genus Eufriesea contains about 60 medium-sized species distributed throughout the Neotropical region. Although no subgeneric classifications exist for this genus, several species' groups were identified by Kimsey (1982) based on subtle morphological characters, including tongue length, integument coloration, genitalic characters and head shape. The genus Euglossa is the most diverse within the tribe, containing approximately 110 species of primarily small bees distributed throughout the Neotropics. This genus is currently subdivided into six subgenera based on external morphological characters (Moure, 1969; Dressler, 1978a, b, 1982a, b, c, d), as well as genitalic characters (Ospina-Torres, Parra & Gonzalez, 2006). Bembé (2007) used morphological characters to infer the phylogenetic relationships among species of the cordata group in the subgenus Euglossa s.s.Darveau et al. (2005) used a single mitochondrial marker to infer the relationships among 32 species of orchid bees, including 24 in the genus Euglossa. Several new species of Euglossa have been described in recent years (Bembé, 2004, 2007; Roubik, 2004; Ramírez, 2005, 2006; Parra, Ospina-Torres & Ramírez, 2006; Rasmussen & Skov, 2006; Hinojosa-Diaz & Engel, 2007; Nemésio, 2007, 2009a, b), and additional taxa probably remain undescribed, especially from poorly sampled areas of South America (e.g. Colombia, Ecuador and Peru). The genus Eulaema contains 15–20 species of large bees (Fig. 1) that are also distributed across the Neotropics. Eulaema is subdivided into two morphologically distinct subgenera (Moure, 1950, 2000). Oliveira (2006) used morphological characters to infer phylogenetic relationships among species, which supported the monophyletic status of the subgenera mentioned above. Finally, Exaerete contains seven cleptoparasitic species that use nests of either Eulaema or Eufriesea as their hosts. Species in this genus are arranged into two species' groups based on morphological traits (Kimsey, 1979; Anjos-Silva, Engel & Andena, 2007).
In summary, the tribe Euglossini encompasses a diverse, but comparatively well-studied group of bees that inhabits the New World tropics. However, despite the widespread interest and extensive literature concerning natural history, behavioural ecology, taxonomy, chemical ecology and population biology – as well as the suitability for comparative study – a comprehensive phylogenetic analysis of the entire tribe is still lacking. Here, we investigate the internal phylogenetic relationships of orchid bees, and explore their morphological evolution, time of divergence, lineage diversification and historical biogeography using comparative methods.
We provide an extensive sampling of the tribe Euglossini, including all five genera, all eight subgenera, 25 of the 26 species' groups and a total of 126 of the approximately 200 species. We sampled nine outgroup taxa, including all corbiculate bee tribes, and two noncorbiculate genera in the family Apidae (Centris and Epicharis). Our DNA matrix contained a total of 146 terminals. Most specimens were collected between 2002 and 2006 from 12 Neotropical countries (Supporting Information Appendix S1). Tissue samples (thoracic muscle or forelegs) or whole specimens were preserved in 100% ethanol and kept at −20 °C for DNA extraction. For each sample, we sequenced a total of ∼4.0 kb from fragments of four different loci, including the mitochondrial protein-coding gene cytochrome oxidase (CO1, 1.2 kb) and the nuclear protein-coding genes elongation factor 1-α (EF1-α, F2 copy, ∼1.2 kb), arginine kinase (ArgK, ∼0.7 kb) and RNA polymerase II (Pol-II, 0.8 kb). Voucher specimens for the sampled taxa (listed in Supporting Information Appendix S1) are currently deposited in the Museum of Comparative Zoology at Harvard University.
We extracted genomic DNA from individual bees using the Qiagen DNA Extraction Kit (Qiagen Inc., Valencia, CA, USA), and ran polymerase chain reactions (PCRs) on a Bio-Rad DNA Engine Dyad Peltier thermal cycler (Bio-Rad Laboratories Inc., Hercules, CA, USA). All reactions were prepared in volumes of 25 µL with 2.5 mmol L−1 MgCl2, 2.5 mmol L−1 PCR buffer and Taq polymerase (Qiagen Inc.) using available primer sets (Danforth et al., 2004a). We purified PCR products by incubating samples at 37 °C for 35 min with the Escherichia coli enzyme exonuclease I (New England Biolabs, Hanover, MD, USA), and subsequently raised the temperature to 80 °C for 20 min. We performed cycle sequencing using a BigDye™ Terminator v3.1 Ready Reaction Cycle Sequencing Kit (Applied Biosystems Inc., Foster City, CA, USA). Samples were sequenced on an Applied Biosystems Inc. 3100 Genetic Analyzer using forward and reverse strands. Complementary strands were assembled with the software package Sequencher v4.2 (Gene Codes Corp., Ann Arbor, MI, USA).
We assembled a DNA matrix containing all four loci using the software packages MacClade v4.06 and Mesquite v2.01 (Maddison & Maddison, 2003). We implemented parsimony analyses in the software package PAUP* v4.10b (Swofford, 2003) by assuming unordered transitions and weighting all characters equally. Heuristic tree searches consisted of 100 random addition sequences, using the tree bisection–reconnection (TBR) swapping algorithm; bootstrap support values were estimated via nonparametric bootstrapping in 1000 replicates. We implemented Bayesian phylogenetic analyses in the software package MRBAYES v3.1.1 (Ronquist & Huelsenbeck, 2003). Tree searches were performed assuming both single and multiple models of sequence evolution for each locus, and Markov chain Monte Carlo (MCMC) searches were made for 10 million generations, sampling every 1000 generations, for a total of 10 000 trees. We estimated model parameters during runs, and estimated Bayesian posterior probabilities as the proportion of trees sampled; the trees obtained in the first one million generations were discarded.
Molecular clock calibration
We performed molecular clock analyses via penalized likelihood (PL) using the software package r8s v1.71 (Sanderson, 2002) and via mean path length (MPL) using the software package PATHd8 v1.0 (Britton et al., 2007). Unlike PL, the MPL method does not assume autocorrelation between parent and child branches, and thus assumes heterogeneous substitution rates across lineages. We calibrated a molecular phylogenetic tree obtained via Bayesian analyses using both fossil and biogeographical data. Two fossil orchid bees are known from Miocene Dominican amber deposits, and several corbiculate bees are known from amber and shale deposits elsewhere. The amber-preserved orchid bee Euglossa moronei (Engel, 1999a) (15–20 Myr old) (Iturralde-Vinent & MacPhee, 1996) bears synapomorphic characters that place it unambiguously within extant lineages of Euglossa. We used its age as a minimum age constraint for the genus Euglossa. The other known fossil orchid bee, Paleoeuglossa mellisiflora (Poinar, 1998), cannot be assigned unequivocally to an extant orchid bee lineage (D. W. Roubik & S. Ramírez, pers. observ.), and thus we did not include it in this analysis. We calibrated our molecular clock analysis with other corbiculate bee fossils, including the giant honey bee Apis lithohermaea from the Pliocene–Pleistocene (14–16 Myr old) olivine basalts of Japan (Engel, 2006), and the stingless bee Proplebeia dominicana from Miocene (15–20 Myr old) lignite, sandy clay beds of Dominican amber deposits and Chiapas amber (Wille & Chandler, 1964; Camargo, Grimaldi & Pedro, 2000). In addition, because orchid bees are restricted to the tropical forests in the New World, we used the Gondwana break-up between tropical America and Africa (∼70–100 Mya; Smith, Smith & Funnell, 1994; Danforth et al., 2004b) as a maximum age constraint for the node connecting Euglossini and Apini. In fact, the fossil record suggests that orchid bees have been restricted to the New World Tropics since the Eocene (Engel, 1999a, b; 2001). On the other hand, extant members of the genus Apis are largely distributed in the Old World, but a recent fossil discovery revealed that honey bees were present in North America (Nevada) during the Miocene (Engel, Hinojosa-Díaz & Rasnitsyn, 2009). We note that molecular clock calibrations based on geographical distributions should be used with caution. The estimated age for corbiculate bees in our analysis is congruent with that suggested by Hines (2008). We estimated branch lengths using the majority rule consensus tree obtained via Bayesian analyses optimized with likelihood using the model of sequence evolution GTR + I + Γ.
Analysis of lineage diversification
We investigated the timing and patterns of diversification in orchid bees with lineage-through-time (LTT) plots. LTT plots depict the cumulative (log) number of nodes as a function of time. We estimated LTT plots on the basis of fossil-calibrated molecular clock chronograms obtained via both PL and MPL. In addition, we tested whether significant changes occurred in the diversification rates throughout the evolutionary history of orchid bees by implementing various likelihood-based statistical methods in the software packages APE v2.0 (Paradis, Claude & Strimmer, 2004) and GEIGER v1.2-02 (Harmon et al., 2008).
Analysis of historical biogeography
Extensive studies of plants and animals have identified recurrent areas of endemism across the Neotropical region. We assigned taxa to biogeographical regions of endemism described by Morrone (2006), and coded the distribution (presence/absence) of orchid bee taxa into discrete geographical regions by surveying both the literature (Kimsey, 1982; Ramírez et al., 2002; Roubik & Hanson, 2004; Nemésio & Silveira, 2007; Nemésio, 2009a, b) and museum collections. For consistency, we used the same nomenclature of biogeographical regions as in Kimsey (1982). In addition, we coded the altitudinal distribution of orchid bee taxa as either predominantly lowland (< 800–900 m), montane (> 900 m) or widespread (Ramírez et al., 2002; Roubik & Hanson, 2004). We analysed the two distributional datasets in a model-based likelihood framework using the software package LAGRANGE v2.0.1 (Ree & Smith, 2008). This flexible platform permits the formulation of dispersal–extinction–cladogenensis (DEC) models by specifying instantaneous transition rates between geographical ranges (or altitudinal zones). Unlike previous methods used in the analysis of historical biogeography (e.g. DIVA), LAGRANGE can incorporate information on past geological events (e.g. land bridge connections, range uplifts) and time (branch lengths), and thus provides a more robust framework for the reconstruction of ancestral areas (Clark et al., 2008).
We built two separate DEC models. In the first model (M0), we allowed all transition probabilities between geographical ranges and altitudinal zones to be equally likely at all times. In the second model (M1), we restricted the geographical ranges and dispersals to reflect geological events that may have affected either vicariance or dispersal events. In M1, lineages were allowed to disperse to and from the Panama endemic zone + Choco region only after the isthmus of Panama was formed gradually as an archipelago at 3.1–3.7 Mya (Coates & Obando, 1996; Coates et al., 2004). Likewise, lineages were allowed to disperse to and from the Andes after the major Andean uplift took place at 2.7 Mya (Gregory-Wodzicki, 2000). The remaining instantaneous transition probabilities were otherwise left unmodified.
Analysis of morphological diversification
To investigate the morphological diversification of functional traits in orchid bees, we used a subset of the species included in the phylogenetic analysis. At least three individuals from a total of 75 orchid bee species (Supporting Information Appendix S2) were measured from specimens deposited in the Museum of Comparative Zoology (Harvard University). We measured total body length, body width (intertegular distance), body mass (dry weight) and tongue length (basal to distal prementum length). We investigated morphological disparity by implementing disparity-through-time (DTT) plots (Harmon et al., 2003). This method calculates dissimilarity as the mean square (pair-wise) Euclidian distance among points in morphospace for subclades relative to the variance of the entire clade. Disparity is calculated for individual nodes by moving up the phylogeny from the root of the tree. This methodology was first implemented by palaeontologists to investigate morphological divergence across geological time scales, and has the advantage of being insensitive to sample size (Foote, 1997). To compare the observed empirical morphological disparity of traits (analysed separately and simultaneously) with respect to a null model, we performed 1000 simulations using both Brownian and speciation models of character evolution. Then, we averaged the simulated values for each node, plotted them against the relative lineage divergence times, and calculated the difference between the empirical and simulated curves by estimating the area between the curves. Brownian models of character evolution assume that a continuous character changes along the branches of a tree following a Brownian motion process, whereas speciation models of character evolution assume that changes on continuous characters occur during branching processes only (Pagel, 1999). These analyses were carried out using the package GEIGER v1.2-02 (Harmon et al., 2008).
Our maximum parsimony and Bayesian inference analyses – based on four loci – resolved most phylogenetic relationships within and between orchid bee lineages. Both methods yielded informative phylogenetic trees (Fig. 3) that were congruent with each other under different optimization schemes and model parameters. Our parsimony analysis yielded 552 shortest trees (TL = 6899), with a strict consensus almost identical to the result of the Bayesian analysis (Fig. 3). The Bayesian phylogenetic analysis also yielded well-supported trees that varied little regardless of whether we used a single (GTR + Γ + I) or partitioned (locus-specific) model of sequence evolution (Supporting Information Table S1). All our combined phylogenetic analyses returned all genera as monophyletic, Exaerete as the sister to the remaining genera, Eufriesea as sister to Euglossa + Eulaema + Aglae, and Aglae clustered within Euglossini (never basal). The relative positions of Eulaema and Eufriesea were inconsistent, depending on the phylogenetic method and model of sequence evolution used.
Our phylogenetic analyses indicated that four of the six subgenera recognized on the basis of morphological characters within Euglossa (Euglossa s.s., Euglossella, Dasystilbe and Glossuropoda) are monophyletic, whereas the remaining two (i.e. Glossura and Glossurella sensuRamírez et al., 2002) constitute paraphyletic groups (Fig. 3). Recently, Moure et al. (2007) and Faria & Melo (2007) proposed the transfer of several species from Glossurella to Glossura, thus making both subgenera more congruent with our phylogenetic results. However, Glossura remains paraphyletic because of the position of Glossuropoda (Fig. 3). In addition, the subgenus Glossurella was not recovered as monophyletic, but this was the result of unresolved nodes (Fig. 3). The remaining species' groups and subgenera in Euglossa, Eulaema and Exaerete were also recovered as monophyletic, but this was not the case for the genus Eufriesea, where most species' groups were recovered as paraphyletic (Figs 3, 4).
Divergence times and diversification patterns
The results from our PL and MPL analyses indicated that Euglossini shared a most recent common ancestor during the Miocene–Eocene periods (27–42 Mya), depending on whether we use the oldest or youngest ages of the fossil calibrations (Supporting Information Table S2). Time estimates obtained via PL and MPL were comparable (Fig. 5).
Our analysis of lineage diversification using LTT plots indicated that orchid bee lineages accumulated rapidly and steadily during the Miocene (15–22 Mya) (Figs 5, 6). The shape of the curve describing lineage accumulation did not change significantly under either PL or MPL, but the slope differed depending on which fossil calibrations were used (Fig. 7). We detected various nodes in which diversification rates increased. Most of these nodes were concentrated in Euglossa.
The probability of detecting changes in the diversification patterns of orchid bees suggested that a model of constant diversification (model A) had a log-likelihood of −366.764 [Akaike information criterion (AIC) = 735.527; δ = 0.183]. A model in which diversification constantly increased or decreased following the Weibull law (model B) had a log-likelihood of −359.083 (AIC = 722.167, α = 0.168, β = 1.30). A model in which diversification abruptly changed with a breakpoint (model C) at 4.0 Mya had a log-likelihood of −390.523 (AIC = 785.045; δ1 = 0.157, δ2 = 0.079). To discern among these alternative hypotheses, we performed likelihood ratio tests, which suggested that diversification following a Weibull law (model B) gave a significantly better fit to our data than the other models (χ2 = 15.361, d.f. = 1, P = 1 × 10−4). Our parameter estimates for model B suggested that β > 1, and thus the diversification of orchid bees has decreased towards the Recent (Supporting Information Table S3). We also carried out a relative cladogenesis test for all slices through the tree, and found that several nodes in the genus Euglossa exhibited an increased diversification rate, even when applying a Bonferroni correction (Fig. 6).
We identified seven areas of endemism for orchid bee taxa across the Neotropical region: (1) Central America; (2) Panama Endemic Zone; (3) Choco; (4) Andes; (5) Amazon; (6) Atlantic Forest; and (7) the Paraguayan Corridor (Fig. 8). The Amazon region exhibited the highest species' diversity (68 species), and the Paraguayan Corridor the lowest (17 species). On average, each species spanned 2.41 ± 1.77 geographical regions. Our dataset contained 56 species endemic to a single region, and eight widespread species whose distribution spanned six or more regions (Fig. 8). Our likelihood analysis of ancestral geographical ranges based on both unrestricted (M0) and restricted (M1) DEC models yielded identical estimates of the global likelihood at the root node (ln L = −438.7). Similarly, dispersal rates (d = 0.05635) and extinction rates (e = 0.01082) did not vary between DEC models M0 and M1. Recurrent range splits between distinct geographical regions were observed. In five cases, single species restricted to Central America were recovered as sister to larger clades that were either restricted to South America or had widespread distributions (Fig. 8; nodes highlighted with red circles). In addition, we observed recurrent range splits connecting the Amazon and the Paraguayan Corridor, the Amazon and the Atlantic Forest, and the Amazon and the Andes (Fig. 8).
Our altitudinal dataset contained 57 lowland species, 23 montane species and 43 widespread species (Fig. 8). The likelihood range reconstruction of altitudinal ranges using DEC models M0 and M1 yielded identical estimates for the global likelihood at the root node (ln L = −152.6), dispersal rates (d = 0.1066) and extinction rates (e = 8.136 × 0−4). Our analysis using the DEC model M1 yielded 14 montane monophyletic lineages (Fig. 8). In all but two cases (Euglossa ioprosopa and E. laevicincta), montane lineages either gave rise to additional montane lineages, or failed to diversify further. On average, montane clades contained two lineages per clade. Montane colonization events coincided with either geographical range shifts from the Amazon to the Andes, or included nonrange shifts within Central America (Fig. 8; grey circles). The species E. ioprosopa and E. laevicincta were unambiguously reconstructed as secondary lowland colonization events derived from exclusively montane lineages (relative probabilities of 0.98 and 0.99, respectively; Fig. 8). On average, the ages of the most recent common ancestor of exclusive montane lineages were 3.88 ± 2.7 Myr and 5.9 ± 4.2 Myr, depending on whether we used the younger or older set of fossil calibrations. Among the 14 montane lineages recovered in our analyses, three were reconstructed as derived from widespread lineages and 10 from exclusive lowland lineages (Fig. 8).
We obtained DTT plots for both individual morphological traits and all traits simultaneously. The results showed substantial variation in the average subclade disparity among traits (Fig. 9). Body width (as measured via intertegular distance) had the lowest average trait disparity with respect to simulated null models of character evolution (inferred using Brownian motion and speciation models of character evolution), as indicated by the lowest values of the area between the observed and null curves (Supporting Information Table S4; Fig. 9). The average empirical disparity calculated for tongue length was lower between subclades (early on during diversification), but showed a marked disparity burst against both Brownian and speciation null curves at the middle of the time scale (Fig. 9). This time period corresponds well with the origin of the species-rich and morphologically diverse genus Euglossa. Moreover, the observed disparity in tongue length decreased towards the tip branches, where the average dissimilarity was substantially lower than predicted by both Brownian and speciation null models (Fig. 9). Conversely, body mass had an intermediate average subclade disparity, as indicated by the close proximity of both the empirical and simulated curves (Fig. 9). However, as with the patterns observed for tongue length, there was a spike in body mass disparity around the middle of the relative time scale. Although the observed disparity in body mass decreased below the null closer to the tips, a large disparity burst was recovered (Fig. 9).
When we used all the morphological variables simultaneously to calculate DTT, the observed average trait disparity accumulated at a slower rate than expected by the Brownian null model, but faster than suggested by the speciation null curve (Fig. 9). Finally, we observed an inverse relationship between both body mass and body size and the species' diversity within each genus (Fig. 10).
The phylogenetic analyses presented here spanned the entire tribe Euglossini, including all genera, all subgenera, most species' groups and about two-thirds of the known species. Our findings shed new light on the phylogenetic relationships, time of divergence, morphological diversification and historical biogeography in orchid bees. Our analyses illustrate the complex interplay between geological events, lineage diversification and trait evolution in the New World tropics.
Orchid bee systematics
Our results indicate that the cleptoparasitic genus Exaerete is sister to the remaining orchid bee genera. Because a reversal from nest-parasitism to nest-building is unlikely (Danforth, 2002; Michener, 2007) and cleptoparasites depend on nest-building species to lay their own eggs, it is unlikely that the most recent common ancestor of Euglossini was a nest-parasite. Thus, our results require either the existence of basal nest-building orchid bee lineage(s) that went extinct, or stem lineage(s) leading to extant Exaerete that switched from nest-building to cleptoparasitism. We speculate that the former scenario is more likely.
Previous systematic studies have addressed the relationships among orchid bee genera, but none of these have suggested the scenario proposed here (Fig. 2). Our results differ from those using morphological characters (Kimsey, 1987; Engel, 1999a; Oliveira, 2006), particularly in the relative position of Exaerete (basal in our study). Such a discrepancy most probably stems from the differences in taxonomic sampling and the number of characters used (< 37). Notably, our results are incongruent with the generic relationships proposed by Michel-Salzat, Cameron & Oliveira (2004). Their study, which was based on both previously available morphological characters and four gene fragments, indicated that Aglae is sister to the remaining orchid bee genera. Although our molecular dataset overlapped partially with two of the gene fragments used by them, none of our gene partitions or combined analyses returned such topology. This incongruence may be explained by their limited sample size (18 species, half in the genus Eulaema).
Our phylogenetic results support previously proposed hypotheses concerning the internal relationships within orchid bee genera, including those for Euglossa, Eulaema and Exaerete. The internal relationships recovered for Euglossa correspond well with several infrageneric groups identified previously on the basis of external morphology (Dressler, 1978a, b). For instance, Euglossa s.s. and all of its species' groups (except E. tridentata) were recovered as monophyletic. Dressler (1978a) defined Euglossa s.s. to include male-specific synapomorphic characters, such as the rhomboid hind tibiae and the comma-shaped anterior tuft of the mid-tibia. Similarly, our results support the monophyletic status of the subgenus Euglossella (+ E. villosa), which was recovered as sister to the rest of Euglossa. The subgenus Euglossella contains some of the most atypical species in the genus, including E. decorata, E. perfulgens and related species, which constitute the only Euglossa s.l. with a brown or black integument that resembles the stingless bees of the genus Melipona (Hinojosa-Diaz & Engel, 2007). Euglossa villosa constitutes the monotypic subgenus Dasystilbe, and is restricted to Central America and Mexico. This species was recovered as sister to the subgenus Euglossella using a locus-specific model of sequence evolution, but as sister to the rest of Euglossa using a single (GTR + I + Γ) model of sequence evolution for all loci. Our phylogenetic analyses also suggest that the subgenera Glossura and Glossurella in the genus Euglossa are paraphyletic. The subgenus Glossura is characterized as containing large bees with long mouthparts, labrum longer than wide and biconvex scutellum (Moure, 1967, 1969), and this was later revised to include species with long, slender bodies and elongated mouthparts (Dressler, 1978a). Glossurella was erected by Dressler (1982a, b) to contain both small and large bees with relatively long tongues, slender body sizes and semicircular depressions on the second sternite of males. Interestingly, both Glossura and Glossurella are not reciprocally monophyletic. Instead, our results suggest that all species from each group possessing a relatively large body size and tongue length exceeding the body length fall into a single, highly supported monophyletic clade (Fig. 3). Moreover, this clade is further supported by the bilobed gonostylus of the male genitalia (Ospina-Torres et al., 2006). The remaining species in the subgenus Glossura, which constitute the bursigera species' group (Dressler, 1982c), were recovered as paraphyletic.
The two subgenera proposed previously for Eulaema were recovered as monophyletic, except for the placement of the species pair composed of E. peruviana and E. speciosa. These results concord well with previously proposed hypotheses based on morphological characters (Moure, 2000; Oliveira, 2006), and are congruent with the intermediate position of E. peruviana between the two subgenera (Oliveira, 2006).
The phylogenetic relationships within the cleptoparasitic Exaerete reflect previous hypotheses about the internal relationships of the genus (Engel, 1999a; Anjos-Silva et al., 2007). Exaerete trochanterica was formerly included in the frontalis species' group (Kimsey, 1979), but our analyses concord with its placement in the dentata group (Anjos-Silva et al., 2007).
Most of the species' groups proposed for the genus Eufriesea were not recovered as monophyletic, which is not surprising given the taxonomic difficulties in this genus and the intergradation of species' groups. In fact, species' groups were originally conceived to aid in species' identifications and not as potential subgenera (Kimsey, 1982). Only a handful of the species' pairs identified by Kimsey (1982) were recovered as closely related in our phylogenetic analyses. In addition, a similar pattern to that found in Euglossa was obtained for Eufriesea. The first branching of the genus corresponds to that between the species' pair Eufriesea caerulescens + E. micheneri (both restricted to Central America and Mexico) and the rest of the genus, which is distributed widely across the Neotropical region.
The molecular clock analyses presented here suggest that orchid bees are a relatively young group, with all extant lineages sharing a most recent common ancestor during the Miocene–Eocene periods, approximately 27–42 Mya. Our age estimates, combined with the species' diversity for each genus, indicate that orchid bees diversified rapidly, particularly at 15–20 Mya.
Our analysis of lineage diversification also indicates that net diversification rates in Euglossini declined through time, as implied by our corrected estimates of the gamma statistic (Paradis, 1997). Numerous recent studies have used species-level molecular phylogenies to explore net diversification rates (speciation rate–extinction rate) and, in most cases, a decrease in the net diversification rate has been observed through time. This pervasive pattern has been attributed to density-dependent diversification processes that may result from a decrease in speciation rates (in the absence of extinction rates) rather than increasing extinction rates (Rabosky & Lovette, 2008). Our study adds weight to this emerging paradigm.
All species of Exaerete are nest-parasites of species of Eulaema or Eufriesea. However, our results suggest that extant lineages in this genus shared a common ancestor before the appearance of these current hosts. Not only has nest-parasitism in Exaerete potentially been stable for 20–25 Myr, but our data also suggest that ancestral host lineages of parasitic orchid bees may have gone extinct. Both Aglae (one species) and Exaerete (seven species) have the lowest species' diversity in the tribe, which may imply that cleptoparasitism imposed constraints on net diversification, as it has appeared to do in other taxa (Wiegmann, Mitter & Farrell, 1993; Pierce et al., 2002). This hypothesis deserves further examination with a more robust phylogenetic sampling of groups across all bees.
Our results also suggest that a negative correlation exists between body size and the number of species per lineage independent of age (Fig. 9). For instance, Aglae (one species), Eulaema (15–20 species) and Exaerete (seven species) have the largest body size among orchid bees, whereas Eufriesea (60 species) and Euglossa (> 100 species) are composed of mostly medium and small bees, respectively. This could be explained by the fact that larger species have greater dispersal capabilities (Janzen, 1971; Dick et al., 2004), and thus could be less prone to geographical isolation, extinction or allopatric speciation. To further explore the role of geography in determining lineage diversification, we mapped the geographical distributions of orchid bees and reconstructed their patterns of historical biogeography using explicit models of dispersal, extinction and cladogenesis.
Historical biogeography and morphological diversification
Our analysis of historical biogeography provides insights into the impact of geological events on the diversification patterns of orchid bees, and reveals a complex interaction between dispersal, local extinction and vicariance at a continental scale. The implementation of explicit DEC models had negligible effects on our estimates of global dispersal and extinction rates, which we interpreted as evidence of congruence between our divergence time estimates and the imposed model constraints (Clark et al., 2008). We observed low local extinction rates (e = 0.01082) relative to dispersal rates (d = 0.05635), which may suggest that extant geographical distributions and diversity patterns in orchid bees were influenced strongly by long-distance dispersal across geographical barriers.
Our analysis of ancestral range reconstruction revealed episodic dispersal events across regions that were disjunct during the Miocene–Pliocene. For example, the most basal node in the genus Eufriesea corresponds to a split between Central American and widespread South American lineages. Similarly, we observed several cases in which early divergences in the genus Euglossa coincided with splits between Central American lineages (including Euglossa jamaicensis; the only described strictly insular species among orchid bees) and South American or widespread lineages (red circles in Fig. 8). Assuming that extant geographical distributions reflect ancestral occurrences, it is tempting to suggest that some early diversification in orchid bees took place after populations moved across the Isthmian archipelago, or even across the Greater and Lesser Antilles, or a proto-Antillean archipelago, with subsequent dispersal and diversification events during the closing of the Panama isthmus at 2.7 Mya (Coates et al., 2004). A similar scenario has been suggested for other groups of animals, including stingless bees (Camargo, Moure & Roubik, 1988). Our likelihood analysis of ancestral range reconstruction also indicates that many lineages of Euglossa diversified within lowland Amazonian forests, and then dispersed across the proto-Antillean Archipelago to give rise to taxa endemic to Mexico. Because we coded the presence of taxa in the Amazon region as binary characters, we ignored distinct distributions within this region. Thus, our dataset is not suited to investigate the fine-grained patterns of historical biogeography within Amazonia (e.g. Dick & Heuertz, 2008; Miller et al., 2008; Santos et al., 2009). Future studies should take into account known areas of endemism for Amazonian orchid bees to investigate the potential impact of Miocene marine incursions (Rossetti, Toledo & Goés, 2005) or Pleistocene forest fragmentation (Mayle, 2004) in the diversification of Euglossini.
We coupled our likelihood estimates of altitudinal ancestral distributions with fossil-calibrated divergence times, and found that montane lineages have diverged from lowland sister clades in the past ∼4–8 Myr. This time range corresponds closely with the estimated time frame of the northern Andean uplift. In fact, approximately 50% of the current Andean altitude was achieved in the past 10 Myr (Gregory-Wodzicki, 2000). Likewise, several Mesoamerican mountain ranges were formed during the Miocene, including the Sierra Madre Occidental and the Neovolcanic Axis (Becerra, 2003 and references therein). An exception to this pattern corresponds to the montane Central American and Mexican species Euglossa villosa, which, based on our clock estimates, diverged approximately 13–17 Mya from its Amazonian sister clade, the subgenus Euglossella. Interestingly, our analyses also indicate that the colonization of montane habitats by orchid bees was relatively rare, and suggest that montane lineages failed to diversify significantly. In contrast with recent studies of vertebrates (Brumfield & Edwards, 2007; Santos et al., 2009), our results indicate that the colonization of the Amazon by Andean montane lineages was exceedingly rare.
Orchid bees exhibit extensive overlap in the functional traits measured across taxa. Our analysis of morphological diversification suggests that changes in body size and tongue length were small between closely related species, and thus were recovered as phylogenetically conserved. Moreover, our analysis suggests that most of the change in body size and tongue length occurred early in the diversification of orchid bees, synchronously with the origin of diverse genera. Although tongue length (Borrell, 2005) and body size (Roubik, 1989, 1992) have been shown to play a role in resource utilization in bees, our analyses suggest that such traits were not involved in the speciation of orchid bees, at least not in the recent past. However, future studies using more fine-grained analyses should consider morphological differentiation in the context of geographical distributions and community assembly.
Euglossine bees account for an important fraction of the pollination services in the Neotropical region (Janzen, 1971; Dressler, 1982a; Williams, 1982; Ramírez, 2009). For instance, several angiosperm plant families, including Costaceae, Euphorbiaceae, Orchidaceae and Zingiberaceae, are known to contain lineages that either rely heavily or depend exclusively on male and female euglossine bees for pollination (Dressler, 1982a). However, the potential impact of such plant–pollinator mutualisms on the evolutionary history of diverse angiosperm groups remains largely unexplored. Orchid bees pose an ideal case to investigate fundamental questions about the role of mutualistic interactions in the origin, evolution and diversification of diverse tropical lineages.
We thank numerous individuals and their respective institutions who provided invaluable assistance in the acquisition of specimens used in this study. We are particularly indebted to T. Arias, R. Ayala, J. Bermudez, S. Cappellari, T. Eltz, C. Garófalo, J. González, A. Link, M. Lehnert, M. Lopez-Uribe, D. McKenna, C. Moreau, J. Neita, L. Packer, A. Parra, H. T. Quental, H. Ramírez, C. Rasmussen, D. Romo, A. Roncancio, T. Samper and C. Yurrita. Funding was provided to SRR and NEP by the National Science Foundation (NSF-DDIG DEB-0608409 and NSF DEB-0447242), and to SRR by the Department of Organismic and Evolutionary Biology (Harvard), the Putnam Fund and Goelet Fund (MCZ, Harvard), the David Rockefeller Center for Latin American Studies (Harvard), the Society of Systematic Biologists Awards for Graduate Student Research, and the Fondo Colombia Biodiversa (FAAE). CS received support from the Danish Research Agency, the Fulbright Commission, the Organization for Tropical Studies (Emily Foster grant), and the University of Florida. Samples were collected (and exported) under permit no. 550 from Colombia (Ministerio del Medio Ambiente), GGVS-309-2003 from Costa Rica (Ministerio de Recursos Naturales, Energía y Minas), AABSX1404311 from Mexico (SEMARNAT) and 005141 from Peru (AG-INRENA). We thank C. Hernandez for assistance with morphometric data. Beth Pringle, Benjamin Bembé, Thomas Eltz and two anonymous reviewers provided useful comments on earlier versions of the manuscript.