Reliable estimates on the ages of the major bee clades are needed to further understand the evolutionary history of bees and their close association with flowering plants. Divergence times have been estimated for a few groups of bees, but no study has yet provided estimates for all major bee lineages. To date the origin of bees and their major clades, we first perform a phylogenetic analysis of bees including representatives from every extant family, subfamily and almost all tribes, using sequence data from seven genes. We then use this phylogeny to place 14 time calibration points based on information from the fossil record for an uncorrelated relaxed clock divergence time analysis taking into account uncertainties in phylogenetic relationships and the fossil record. We explore the effect of placing a hard upper age bound near the root of the tree and the effect of different topologies on our divergence time estimates. We estimate that crown bees originated approximately 123 Ma (million years ago) (113–132 Ma), concurrently with the origin or diversification of the eudicots, a group comprising 75 per cent of angiosperm species. All of the major bee clades are estimated to have originated during the Middle to Late Cretaceous, which is when angiosperms became the dominant group of land plants.
Bees are among the most important pollinators of flowering plants (angiosperms), of which 78–94% of species rely on animal pollinators . Angiosperms include 352 000 species  and dominate terrestrial ecosystems. Perplexed by what appeared to be the sudden appearance of modern flowers in the fossil record, Darwin referred to this great diversity as an ‘abominable mystery’ . Angiosperms became dominant in species numbers during the lower Late Cretaceous (especially from the Albian to Turonian age) . The matter of angiosperm success has received much attention and continues to be a highly debated and controversial area of research (summarized in ). It has often been hypothesized that bees arose concurrently with the diversification of flowering plants [5–9] and could have caused rapid and extensive speciation in angiosperms through ethological and floral reproductive isolation [10,11], or by increasing specialization in traits involved in biotic pollination within a given area .
In order for bees to have contributed to the diversity of flowering plants seen today, they must have been present and diversifying at the same time as angiosperms were diversifying. This requires reliable estimates on the dates of divergence events in both groups. Much research has been focused on dating the origins of the major angiosperm clades (reviewed in ). Recent Bayesian analyses using relaxed molecular clock methods provide a wide range of dates for crown group angiosperms ( 221.5–275 Ma (million years ago);  147–183 Ma;  167–199 Ma). Regardless of exactly when angiosperms originated, information from the fossil record clearly shows that a major radiation of angiosperms occurred during the Mid-Cretaceous .
In contrast to angiosperms, studies putting precise dates on the antiquity of the major clades of bees using information from both DNA sequence and fossil data have been lacking. To date, there have only been age estimates using such methods for particular groups of bees [17–32]. Estimation of divergence events in bees has been hampered by a lack of well-supported phylogenies and a scarcity of bee fossils. During the past two decades, great advances in both of these areas have been made, making this an opportune time to incorporate information from phylogenetic hypotheses, the fossil record and rates of molecular evolution, into estimates of the ages of major bee clades. Phylogenetic analyses using molecular data or combining morphological and molecular data have elucidated relationships within and among many bee lineages (reviewed in ). In addition, numerous bee fossils have now been recovered and described [34–45].
Despite the recent publication of new bee fossil taxa, the bee fossil record remains sparse and highly biased towards resin-collecting bees, which represent 61 per cent of the described species of bee fossils  but comprise only 29 per cent of extant bee species. The oldest described fossil bee, Melittosphex burmensis, is from Burmese amber estimated to be 100–110 Myr (million years old) [39,42]. This bee appears to be morphologically intermediate between bees and crabronid wasps (the putative sister group to bees), and does not seem to belong to any extant family of bees. Ohl & Engel  also argue that it may in fact not be a bee, but a predatory wasp. Therefore, this fossil cannot be used as a minimum age estimate for the crown group of bees, as it may represent a stem lineage. The oldest fossil that can be attributed to the crown group of bees with confidence is Cretotrigona prisca from the Late Cretaceous amber of New Jersey . The exact age and phylogenetic affinity of this fossil has been controversial (see calibration node 5 in the electronic supplementary material). However, this fossil is a stingless bee belonging to the crown group or stem lineage of Meliponini [47,48], a tribe within Apidae and not near the base of bee phylogeny. This fossil is probably of Maastrichtian age (approx. 65 Myr) . The lack of crown group bee fossils from the Cretaceous favours the hypothesis of a Tertiary radiation of bees. However, the existence of an approximately 65 Myr stingless bee suggests the possibility that many of the major bee clades were already present by the end of the Cretaceous and that bees radiated during the Cretaceous, coincident with the diversification of angiosperms.
An incomplete fossil record inherently underestimates the ages of clades. One way to estimate divergence events in the absence of a complete fossil record is to use model-based methods that allow for the combination of substitution rate data (obtained from DNA sequence data) and calibration points (obtained from the fossil record). A number of methods have been developed to estimate divergence dates under models that relax the strict molecular clock constraint [49,50]. More recently, a Bayesian Markov chain Monte Carlo (MCMC) method has been introduced for performing relaxed phylogenetics in which the phylogeny and the divergence dates are co-estimated under a relaxed molecular clock using probabilistic calibration priors .
In this paper, we expand the molecular dataset of Danforth et al. , by adding two new protein coding genes and adding 75 more taxa (for a total of seven genes and 168 species). The addition of these new taxa allowed us to use information from the fossil record to time calibrate 15 nodes in the phylogeny and estimate divergence times of the major bee clades using the relaxed phylogenetics approach implemented in Beast . We explore the effect of applying a hard upper bound near the root of the tree on our divergence time estimates and investigate the effect of various tree topologies obtained using different Bayesian methods on our age estimates. We then compare our age estimate for bees with those of angiosperms to see whether there is a correspondence in timing of major divergence events in bees and angiosperms.
(a) Taxon sampling
We sampled representatives from all subfamilies of bees and 91 per cent of tribes (unsampled tribes: Disoglottini, Protandrenini, Nolanomelissini, Protomeliturgini, Perditini) following . Our dataset represents all major lineages of bees with 152 species included. Care was taken to select representatives that would allow us to more accurately place fossil calibration points. As outgroups, representatives from all three subfamilies of Sphecidae and four of the five subfamilies of Crabronidae were included, because bees appear to have arisen from within the ‘spheciform’ wasps  and more specifically to be sister to the aculeate wasp family Crabronidae [55,56] or to have evolved from within Crabronidae . All of the species included in the study along with taxonomic and voucher information are listed in the electronic supplementary material, table S1.
The dataset consists of sequences from two nuclear ribosomal genes (18S, 28S), and five nuclear protein-coding genes (wingless, pol II, opsin, Nak, Ef1α). Previously published sequences were downloaded from Genbank (see the electronic supplementary material, table S1). All new sequences were obtained following standard PCR and sequencing protocols . Primer pairs and PCR conditions for all genes are listed in electronic supplementary material, table S2. All genes were separately aligned in the Lasergene DNAStar software package using ClustalW. Alignments for 28S and 18S were subsequently adjusted by referring to the secondary structure of these genes proposed for Apis mellifera , and those regions that could not be aligned with confidence were excluded from the analysis. Intron regions of opsin and EF1α also could not be aligned unambiguously and were therefore excluded from the phylogenetic analysis. The final number of aligned nucleotides used in the analyses was 5513 (18S: 740 bp; 28S: 881 bp; wingless: 403 bp; pol II: 841 bp; opsin: 456 bp; Nak: 1434 bp; EF1α: 758 bp), and in the final data matix only 3 per cent of cells were missing data.
(i) Phylogenetic analyses
The dataset was divided into four partitions. The two ribosomal genes were placed together into one partition, and the protein-coding genes were combined together and partitioned by codon position (|28S,18S|pos1|pos2|pos3|). The general time reversible (GTR) model , with a proportion of invariable sites (I), and rate variation among sites with four rate categories (G) was assigned to all partitions based on the results of model tests done in jModelTest v. 0.1.1 . Eight independent analyses with four chains each were run using MCMC methods in MrBayes v. 3.1.2 . The number of generations for each run varied from 9 741 000 generations to 17 824 000 generations. The parameter trace files of each run were observed in Tracer v. 1.4.1  to verify that the runs had converged on the stationary distribution, and to decide on the appropriate number of generations to discard as burn-in (individual runs summarized in electronic supplementary material, table S3).
The tree files with the burn-in removed from each run were combined giving a total of 81 396 post-burn-in trees. A maximum clade credibility tree was constructed from 6434 trees subsampled evenly from the combined tree file in TreeAnnotator v. 1.4.8. We chose to do numerous shorter runs instead of a few longer runs because of run time restrictions on the computer clusters at Cornell University's Computational Biology Service Unit http://cbsuapps.tc.cornell.edu/. Also, doing a large number of independent runs from different starting points allowed us to more fully explore tree space.
(ii) Estimating divergence times
We used a Bayesian uncorrelated relaxed-clock (UCLN) model  with multiple calibration points to estimate divergence times in the program Beast v. 1.6.1 . We partitioned the dataset and applied a GTR + I + G model to each partition, as in the phylogenetic analysis described above, and allowed the tree topology to be estimated to accommodate for phylogenetic uncertainty. Under the UCLN model, branch substitution rates are drawn independently from an underlying log-normal distribution. This allows for the rate of evolution to vary among the branches of the tree with no a priori correlation between a lineage's rate and that of its ancestor. The Yule tree prior was used, which assumes a constant per lineage selection rate as recommended for species-level phylogenies . We randomly selected a starting tree from the posterior distribution of trees from the MrBayes analysis and scaled this tree so that it was consistent with all of our calibration points.
The tree was time-calibrated by applying a prior probability on the ages of 15 internal nodes (see the electronic supplementary material, figure S1, table S5). Age estimates were based on paleontological evidence as described in the electronic supplementary material. There are 184 described bee fossils , but we only used fossil taxa which could be confidently assigned to taxonomic groups and nodes represented in our dataset, and for which precise stratigraphic information was available. Uncertainty in the age of the calibration points was incorporated into the analysis by assuming that the prior probability of the node being a certain age follows a lognormal distribution with a rigid minimum bound. Applying a lognormal distribution to our age estimates allows us to assume that the actual divergence event took place sometime prior to the earliest appearance of fossil evidence, but that the age of the node is more likely to be close to the age of the oldest known fossil and less likely to be significantly older. By how much the appearance of a clade predates the age of the first fossil is always unclear. We therefore selected a mean and standard deviation for each calibration point so that the 95 per cent prior probability for the age of the node ranged from two to 25 million years older than the earliest fossil evidence for the group. In a few cases, we were not confident that the fossil was part of the crown group and not the stem group. In these cases, we applied a normal distribution to the calibrating node so that it was just as probable that the node was 10 million years older or younger than the fossil.
Our initial divergence time analysis described above (referred to as analysis 1 below) does not use any hard maximum bounds on the age of any nodes. Unfortunately, this means that there are no hard upper constraints to divergence times which can lead to arbitrarily old age estimates (especially near the root of the tree) while the substitution rate takes on low values . Even when each node has a prior age distribution associated with it, as is the case here, a more explicit prior on the age of root is desirable . We therefore also ran analyses in which we placed a normally distributed prior on the age of the root node of the tree (which represents the divergence between Sphecidae and Crabronidae + bees). The oldest apoid fossils are from the Early Cretaceous and belong to the extinct stem group lineages referred to as Angarosphecidae . The oldest fossils that can be assigned with certainty to extant apoid families come from Myanmar amber . Fossils of Ampulicidae and Crabronidae (Pemphredoninae) obtained from this amber have been assigned to the Upper Albian (approx. 100–110 Myr) of the Early Cretaceous . Taking into consideration the entire apoid fossil record, Grimaldi & Engel  estimate that Sphecidae diverged from other Apoidea during the Early Cretaceous approximately 140 Ma. We therefore ran analyses with the age of the root node sampled from a normal distribution with a mean of 140 and standard deviation of five (referred to as analysis 2 below) to allow for some uncertainty in the age of the root node. Because the tree topology obtained from our Beast analyses differed from that obtained in the MrBayes analysis, we also ran analyses in Beast in which the tree was constrained to have the same relationships as the tree obtained by MrBayes at the subfamily and family levels (referred to as analysis 3 below).
For each analysis, over 200 million post-burn-in generations were accumulated from numerous individual runs which each ranged from 16 million to 100 million generations (individual runs summarized in electronic supplementary material, table S4). The parameter trace files of each run were observed in Tracer v. 1.4.1  to verify that the runs had converged on the stationary distribution, to decide on the appropriate number of generations to discard as burn-in, and verify that the estimated sample size of each parameter was above 200. A maximum clade credibility tree was constructed as described above.
Our MrBayes analysis of the concatenated seven gene dataset results in a well-resolved phylogeny with most nodes having high support (see the electronic supplementary material, figure S1). The monophyly of bees is supported with posterior probability (PP) of 1.0. The melittid subfamily Dasypodainae is sister to all other bees (PP = 0.86). Meganomiinae and Melittinae form a monophyletic group (PP = 1.0) sister to all non-melittid bees (PP = 0.97). The long-tongued bees are monophyletic (PP = 1.0) as are the two long-tongued bee families, Apidae (PP = 1.0) and Megachilidae (PP = 1.0). The short-tongued bee families Andrenidae, Halictidae, Colletidae and Stenotritidae form a monophyletic clade (PP = 0.99) sister to the long-tongued bees. The monophyly of these four families are each recovered with PP = 1.0. Andrenidae is found to be sister to the remaining members of this group (PP = 1.0), and Colletidae and Stenotritidae form a monophyletic group (PP = 1.0). Our study recovers the same phylogenetic relationships among bee families and subfamilies as Danforth et al. .
The maximum clade credibility tree resulting from our initial Bayesian uncorrelated relaxed-clock model with multiple calibration points divergence time analysis (analysis 1) (see the electronic supplementary material, figure S2) showed different relationships than those obtained from our MrBayes analysis (see the electronic supplementary material, figure S1). Our analysis 1 tree showed the initial divergence in bees to be between the short-tongued and long-tongued bees. To see whether this change in topology was because of the calibration points used, an additional divergence time analysis was run in which all internal calibration points were removed and instead a prior on the root node age was placed with a mean ± s.d. of 140 ± 5 Ma. Removal of the internal calibration points did not alter the tree topology from that obtained in analysis 1. To see whether the change in topology seen in analysis 1 was because of uneven taxon sampling from each of the major bee clades, an additional divergence time analysis was run in which taxa were removed so that all major clades had proportional representation to the number of species described in that clade. This modification did not succeed in producing the same topology as that found in the MrBayes analysis. Removing these taxa also did not noticeably alter the estimated age of the common ancestor of extent bees from that estimated in analysis 1. We also ran an analysis in Beast applying a strict clock instead of a relaxed uncorrelated clock and found this tree to differ from both the MrBayes tree and the tree obtained in analysis 1.
When a prior on the root node age was placed in addition to the internal calibration points (analysis 2), extant bees are estimated to be 111 to 130 Myr (see the electronic supplementary material, figure S3). Constraining the tree topology to be congruent with the MrBayes tree (analysis 3) (figure 1) barely altered our estimate on the age of extant bees (113–132 Myr). The estimated age for all major bee clades under the three different analyses are listed in electronic supplementary material, table S6, along with estimates from previous studies.
(a) Effect of root node prior and tree topology on divergence time estimates
Our estimate for the antiquity of bees is much older when no prior is placed on the age of the root node (figure 2). As we move closer to the tips of the tree, placing a prior on the root node has less of an effect on estimated ages and we see an overlap between the 95 per cent HPDs of all three analyses for the ages of the major bee clades. Because there were no hard upper constraints to divergence times in analysis 1, we believe that this led to arbitrarily old age estimates (especially near the root of the tree). Ho & Phillips  pointed out this artefact of estimating old divergence times when no hard upper age constraints are used. Our estimates for the antiquity of bees based on analyses 2 and 3 are also more consistent with previous hypotheses on the age of bees based on extrapolating from the fossil record alone . Our age estimates for all of the major bee clades did not meaningfully differ between analyses 2 and 3, despite the differences in topology. For this reason, we restrict our discussion below to only the analysis 3 age estimates, in which a root node prior is applied and the tree topology is constrained to be congruent with that recovered from the MrBayes analysis.
(b) Estimated ages in comparison to previous clade-specific studies
We estimate that the common ancestor of all extant bees originated approximately 123 Ma (113–132 Ma) and that all extant families were most probably present before the Cretaceous/Paleogene transition (formally referred to as the K/T boundary) and period of mass extinction (except for Stenotritidae) (figures 1 and 2). Our estimated ages based on analysis 3 for the various groups are in large part congruent or have overlapping ranges with previous studies (see the electronic supplementary material, table S6). Exceptions include our ages estimated for various megachilid groups when compared with those estimated by Litman et al. . Because they did not use any strong upper bounded age prior near the root of their tree, their ages are consistent with our ages estimated in analysis 1. Our age for bumblebees is also rather young compared with those estimated by Hines  and Ramirez et al. . The estimated age ranges however overlapped in most cases and our younger dates may be because of incomplete taxon sampling. As more species are sampled within a clade, the branch lengths within that clade become longer and age estimates can consequently also become older. Therefore, our study may be underestimating ages of clades at a lower taxonomic level when compared with clade-specific analyses that include a much higher proportion of species. However, a large-scale study such as ours has the advantage of being able to more fully use information from the fossil record to place multiple calibration points throughout the tree.
Likewise, our much younger age estimate for crown group Meliponini compared with that obtained by Rasmussen & Cameron  may be due in part to our incomplete taxon sampling, but it is more probably because of the difference in how we used the extinct fossil bee Cretotrigona prisca as a calibration point. Rasmussen & Cameron  used this fossil to place a minimum age constraint of 65 Ma on crown group stingless bees. We are, however, not confident that this fossil belongs to crown group Meliponini and therefore made it a priori equally probable that crown group Meliponini is older or younger than 65 Myr. Despite this, the analysis strongly favoured a younger age of 49–59 Ma, indicating that placing a minimum age 65 Ma on this node would be greatly inconsistent with information from the other calibration points. Near et al.  suggested removing inconsistent calibration points from molecular clock analyses, but we prefer not to do so given the sparse nature of the bee fossil record. Instead, we interpret the different age estimates obtained by our study and the Rasmussen & Cameron  study as support for the hypothesis that Cretotrigona prisca is a stem lineage of Meliponini.
(c) Estimated ages in comparison to bee traces in the fossil record
We did not use any paleoichnological data as calibration points in our divergence time analyses even though for some groups of bees, traces of nests and leaf damage thought to have been done by particular bees are from an earlier time period than any actual body fossils. For example, fossil nests from the Late Cretaceous have been attributed to halictid bees  whereas the oldest halictid bee body fossil Halictus savenyei from the Quilchena deposit (Canada) is 53 Myr . Despite not using any ichnotaxa as calibration points, our estimate of 75–96 Myr for Halictidae supports the view that these nest traces could belong to halictid bees. Conversely, leaf fossils from the Mid-Eocene with damage that has been attributed to the work of bees from the tribe Megachilini [69,70] are much older than our estimated age for crown group Megachilini. However, our node representing Megachilini may not be a representative of the common ancestor of all extant Megachilini because in a recent molecular phylogenetic analysis of Megachilidae , Coelioxys and Radoszkowskiana arose from within Megachile. Therefore, our estimated age for crown group Megachilini is probably somewhere along the branch leading to the node uniting Coelioxys and Megachile, which is quite long, and would accommodate a Mid-Eocene or older date for the origin of Megachilini.
(d) Estimated ages in relation to angiosperm divergence events
Our estimate that the common ancestor of all extant bees originated approximately 123 Ma (113–132 Ma) is more recent than the recently estimated ages for the origin of angiosperms based on relaxed fossil calibrated molecular clocks. Smith et al.  estimated that crown group angiosperms originated 182–270 Ma, Bell et al.  141–199 Ma and Magallon  221–275 Ma. Estimates for the age of crown group angiosperm based on the fossil record alone are consistent with a Valanginian origin . The oldest angiosperm macrofossil, Archaefructus lianingensis from the Yixian lake beds of northeast China [71,72] is of Barremian or Aptian age . Some pre-Cretaceous angiosperm-like fossils have been described [74–76], but their assignment to angiosperms remain tentative at best .
Our estimated age for crown group bees, is however, coincidental with the first appearance of tricolpate pollen grains (approx. 125 Ma) which is often equated with the origin of eudicots . Because tricolpate pollen is a unique synapomorphy, is conspicuous in fossil state and has a dense fossil record, it is thought that if eudicots were present prior to the Aptian, there would probably be evidence of tricolpate pollen earlier in the fossil record . Divergence time analyses using relaxed uncorrelated molecular clocks differ somewhat in their estimates on the age of eudicots. Bell et al.  estimate that eudicots originated 123–139 Ma whereas Smith et al.  estimate 128–172 Ma. Magallon  used the age of tricolpate pollen to place a maximum age constraint of 125 Ma on the eudicots, indicating her belief that eudicots are approximately 125 Myr. As noted by Smith et al. , if these older eudicot age estimates are correct, then the appearance of tricolpate grains may instead signal the rise in abundance and geographical expansion of eudicots. Eudicots represent approximately 75 per cent of the diversity of flowering plants , and a large proportion of families are highly dependent on bees for pollination services. Similarly, bees rely heavily on eudicots for pollen and nectar [78,79] and floral oils . The oldest inferred gain of oil production in flowers was in Malpighiaceae (75–64 Ma) . Bees of Centridini and Tetrapedia forage for oil on various Malpighiaceae and we estimate that these bees evolved around that time period (Centridini: 52–87 Ma, Tetrapedia stem lineage: 66–92 Ma).
Our estimated age for bees is also concurrent with an increase in specialized pollination modes within angiosperms. Hu et al.  found evidence of adaptations to permit pollen clumping and more specialized pollination modes during the Mid-Cretaceous angiosperm diversification. They hypothesized that the increase in specialized pollination modes seen during this period was linked to bee pollination. Vamosi & Vamosi  hypothesized that the evolution of specialized pollination modes may have had an effect on angiosperm diversification based on their finding that traits associated with biotic pollination had the greatest impact on diversification in angiosperm families once area was taken into account. Prior to the Mid-Cretaceous and our estimated age for the origin of crown group bees, fossil flowers show characteristics consistent with more generalized forms of pollination by Coleoptera, Diptera and Lepidoptera .
More specific evidence of the dependence of angiosperms on bees dates back to the Turonian. Fossil Clusiaceae from the Late Cretaceous (Turonian) have highly distinctive morphological features that strongly suggest dioecy and resin production . Extant species with resin-producing flowers are often pollinated by bees that collect resin mostly for nest construction and defensive antimicrobial compounds, such as stingless and orchid bees [83–85]. Based on our fossil-calibrated age estimates for the corbiculates, we estimate that stingless and orchid bees were not present during the Turonian, but that extinct corbiculate stem lineages were. Consistent with this hypothesis, extant honey, stingless and orchid bees all collect resin, making it highly probable that the common ancestor of corbiculate bees (and therefore extinct stem lineages) also collected resin.
Further evidence of angiosperms depending on bees for pollination during the Turonian comes from fossils of Ericaceae (Ericales) found in Turonian deposits. These fossil flowers show a complex of characters (e.g. elongate sepals, nectary disc, carpels with separate stigmas, a sympetalous corolla, two whorls of stamens with stamen awns) that are now associated with highly specific, usually apid bee pollinators . We estimate that apids were present during the Turonian and therefore provide further support for the pollination of these extinct Ericaceae by apid bees.
We estimate that crown group bees originated approximately 123 Ma (113–132 Ma), concurrently with the origin or the rise in abundance and geographical expansion of eudicots, which are heavily dependent on bees for pollination. Most of the bee families originated during Mid- to Late Cretaceous when angiosperms became dominant in species numbers. Therefore, our results support the view that bees were present early in the evolution of angiosperm plants and that they could have played an important role in diversification of the eudicots, a group including over 75 per cent of extant angiosperm species.
This work was in part supported by NSF grants DEB 0709956 to S.C. and DEB 0814544 to B.N.D. We thank Andrew Debevec for help with data collection.
- Received November 16, 2012.
- Accepted January 4, 2013.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.