Molecular phylogeny and historical biogeography of the large carpenter bees, genus Xylocopa (Hymenoptera: Apidae)

The biogeographical history of major groups of bees with worldwide distributions have often been explained through hypotheses based on Gondwanan vicariance or long distance dispersal events, but until recently these hypotheses have been very difﬁcult, if not impossible, to distinguish. New fossil data, comprehensive information on Mesozoic and Cenozoic coastline positions and the availability of phylogenetically informative DNA markers now makes it fea-sible to test these hypotheses for some groups of bees. This paper presents historical biogeographical analyses of the genus Xylocopa Latreille, based on phylogenetic analyses of species belonging to 22 subgenera using molecular data from two nuclear genes, elongation factor-1 α ( EF-1 α ) and phosphoenolpyruvate carboxykinase ( PEPCK ), combined with previously published morphological and mitochondrial data sets. Phylogenetic analyses based on parsimony and likelihood approaches resulted in several groups of subgenera supported by high bootstrap values ( > 85%): an American group with the Oriental/Palaearctic subgenera Nyctomelitta and Proxylocopa as sister taxa; a geograph-ically diverse group ( Xylocopa s .l ); and a group consisting of African and Oriental subgenera. The relationships among these three clades and the subgenus Perixylocopa remained unresolved. The Oriental subgenus Biluna was found to be the sister group of all other carpenter bee subgenera included in this study. Using a relaxed molecular clock calibrated using fossil carpenter bees, we show that the major splits in the carpenter bee phylogeny occurred well after the ﬁnal breakup of Gondwanaland (the separation of South America and Africa, 100 Mya), but before important Miocene fusion events. Ancestral area analysis showed that the genus Xylocopa most likely had an Oriental-Palaearctic origin and that the present world distribution of Xylocopa subgenera resulted mainly from independent dispersal events. The inﬂuence of Pleistocene glaciations on carpenter bee distributions is also discussed. © 2002 The Linnean Society of London, Biological Journal of the Linnean Society , 2002, 77 , 249–266.


INTRODUCTION
In a paper published in 1979, Charles Michener comprehensively discussed the origins and evolution of the major bee groups of the world and this review has continued to provide the theoretical framework for interpreting historical biogeographical patterns in bees.A central tenet of this review was that several closely related bee groups (tribes or genera) have disjunct distributions that are best explained by vicariance associated with the breakup of Gondwanaland.A good example is a group of bees of the tribe Paracolletini that only occur in Australia, South Africa and temperate South America (Michener, 1979: 279, 301-302).Other groups of bees, because of their shared occurrence on adjacent continents (for example, related species in the genera Andrena and Halictus in North America and Europe), suggest relatively recent dispersal events.At the time of Michener's review the standard approach to investigate historical biogeography was to use distribution patterns, species richness of areas and relationships between bee groups based on morphology.However, there were few cladistically based phylogenetic trees available for groups with worldwide distributions, making it difficult to interpret their historical biogeographical patterns and distinguish between alternative hypotheses, such as Gondwanan vicariance or intercontinental dispersal events.With the application of DNA data to resolve phylogenies, which have the potential for divergence dates to be estimated using molecular clock techniques, and the availability of comprehensive information on the positions of Mesozoic and Cenozoic coastlines (Smith, Smith & Funell, 1994), it is now possible to re-examine these hypotheses.In the following paper we carry out such an analysis of the historical biogeography of the large carpenter bees, genus Xylocopa Latreille 1802.
Xylocopa consists of about 469 species (Michener, 2000) in 31-51 subgenera, depending on the classification that is followed.In a taxonomic revision of the tribe Xylocopini (Hurd & Moure, 1963) three different genera were recognized: Lestis Lepeletier & Serville, Proxylocopa Hedicke, and Xylocopa Latreille comprising 51 subgenera.In a cladistic analysis of the genera and subgenera based on parsimony analyses of morphological characters, Minckley (1998) concluded that the genera Lestis and Proxylocopa should be considered as subgenera of Xylocopa to avoid paraphyly of this genus, a viewpoint which is adopted in the present paper.Minckley (1998) synonymised several subgenera and erected one new subgenus, recognizing 33 different subgenera.Although his taxonomic decisions were based on the observations that the synonymised taxa grouped together in all phylogenetic analyses, and were supported by one or more synapomorphies, some of the new combinations (e.g.Koptor- tosoma and Mesotrichia ) now link subgenera that are from disjunct geographical regions.To be conservative, and to maintain the geographical resolution among the subgenera, for this paper we will follow the classification of Hurd & Moure (1963).
Xylocopa is distributed over most continents, predominantly in tropical and subtropical climates and occasionally in temperate areas (Hurd & Moure, 1963).Most of the subgenera have restricted distributions that do not cross boundaries of the world's main zoogeographical regions.Some subgenera are endemic to islands: Hoplitocopa (Lesser Sunda Islands, Indonesia), Lieftinckiella (Sulawesi) and Prosopoxylocopa (Madagascar).A few subgenera are distributed over more than one zoogeographical region.For example, most species in the subgenera Neoxylocopa , Schoen- herria and Stenoxylocopa are neotropical, but only a few species, Neoxylocopa mexicanorum Cockerell and N .varipuncta Patton (two out of 50 species), Schoen- herria loripes Smith and S .micans Lepeletier (two out of 33 species), Stenoxylocopa m. micheneri (Hurd) (one out of five species), are found in nearctic regions (Hurd, 1978), and are likely to have originated by dispersal from South America.The subgenus Koptorto- soma , currently with about 130 recognized species, has two centres of diversity, one in Africa, south of the Sahara (Eardley, 1983) and one Oriental (Maa, 1938;van der Vecht, 1953;Lieftinck, 1955Lieftinck, , 1956)), that extends into the Australo-Papuan region (Lieftinck, 1957;Leys, 2000).The same kind of disjunct distribution is found in the subgenus Ctenoxylocopa , whose species are found from southern Africa through the Middle East to Indonesia (Maa, 1970).
The fact that the majority of the Xylocopa subgenera are confined to single zoogeographical regions (Hurd & Moure, 1963), has led to several intuitive explanations about the origin of the genus and two hypotheses have been formulated to explain the present distribution of subgenera (Hurd & Moure, 1963;Hurd, 1978;Michener, 1979).The hypotheses mainly differ in the timing of the origin of the genus.A Gondwanan origin of the genus has been used to hypothesize that distinct subgeneric groups, for example from South America and Africa, are the result of past vicariance events (Gondwana break-up).A post-Gondwanan origin of the genus, however, has been suggested based on observations that pairs of presumably related taxa occur on both sides of the North Atlantic and the Pacific.These taxa are likely to have originated following intercontinental dispersal of ancestral taxa via island hopping or across land bridges (Hurd & Moure, 1963;Michener, 1979).Although carpenter bees are medium to large bees that are capable of flying relatively large distances over water (a species colonized Krakatau shortly after the eruption of this volcanic island; Docters van Leeuwen, 1936), there is no evidence for long distance (intercontinental) dispersal.The few Xylocopa species that occur on pacific islands (Hawaiian islands, Galapagos Islands) most likely arrived there relatively recently through human mediation (Hurd & Moure, 1963;Hurd, 1978).To test alternative biogeographical hypotheses, and to place these hypotheses in an appropriate time frame, e.g.whether disjunct distributions of pairs of related species (or species groups) have evolved through dispersal or vicariance, it is crucial to work with a robust phylogeny (Avise, 1994;Voelker, 1999).
Recently, phylogenetic analyses were undertaken to resolve evolutionary relationships between subgenera of Xylocopa using a molecular data set comprised of sequences of two mitochondrial DNA genes (Leys, Cooper & Schwarz, 2000).Although both studies resolved some major groups of subgenera, the relationships between those groups and other subgenera were not resolved.The main reasons for this lack of resolution were: (1) the lack of phylogenetically informative characters and large number of homoplasious characters in the morphological data set (Minckley, 1998), and (2) an extreme A-T nucleotide bias, leading to homoplasies, particularly at transitions of third codon positions in the mtDNA data set, and the fact that basal divergences in the phylogeny were only supported by a few synapomorphies (Leys et al ., 2000).
We therefore aim to further resolve the phylogenetic relationships among subgenera within the genus Xylo- copa by using sequence data from the two nuclear genes elongation factor-1 α ( EF-1 α ) and phosphoe- nolpyruvate carboxykinase ( PEPCK ).These nuclear genes are known to evolve more slowly than mitochondrial genes and may be useful to resolve deeper divergences in the phylogeny (Friedlander et al ., 1996).In particular, the PEPCK gene has been successfully used to resolve Mesozoic-age divergences within Lepidoptera (Friedlander et al ., 1996), and is therefore expected to be useful for subgeneric analyses of Xylo- copa as well, given the potentially Late Cretaceous/ Early Tertiary origin of Xylocopa (Engel, 2001a).Evidence from Southern analyses, using a Drosophila probe, suggest that PEPCK is a single copy gene in Lepidoptera and Diptera (Friedlander, Regier & Mitter, 1992;Friedlander et al ., 1996).The PEPCK gene includes several introns that might be informative for relationships among recently diverged taxa.EF-1 α has proven useful for phylogenetic analyses of a variety of other insect groups, like swallowtail butterflies (Reed & Sperling, 1999) and bees in the family Halictidae (Danforth, Sauquet & Packer, 1999).Two copies of EF-1 α (F1 and F2) have been found in honeybees ( Apis ) and flies ( Drosophila ) which differ in sequence by approximately 19% and 25%, respectively, and they also differ in the number and location of introns (Danforth & Ji, 1998).
The use of molecular data for inferring phylogenetic relationships opens the possibility to estimate divergence times of nodes when fossil data are available.Fossil carpenter bees likely to be related to the subgenera Xylocopa s. s. , Copoxyla and Nyctomelitta are known from Middle Oligocene and Miocene deposits (Hurd & Moure, 1963;Zeuner & Manning, 1976;Zhang, 1990), and from the Middle Eocene (Engel, 2001a).The use of molecular clocks is often criticised because of accumulating data that rates of molecular evolution may not be constant (Wu & Li, 1985;Muse & Weir, 1992;Sorhannus & Van Bell, 1999).However, recent methods that incorporate models of change of molecular rates in and among branches (Sanderson, 1997;Thorne, Kishino & Painter, 1998;Huelsenbeck, Larget & Swofford, 2000;Kishino, Thorne & Bruno, 2001) make it possible to linearize phylogenetic trees, and estimate divergence times when the age of one or more nodes can be calibrated using fossil data.
The aim of this paper is to combine phylogenetic data, the fossil record and data on positions of palaeogeographical coastlines to test vicariance-dispersal hypotheses, and to synthesize a plausible hypothesis of the historical biogeography of large carpenter bees.

T AXA EXAMINED
In order to study the phylogenetic relationships among the subgenera of the genus Xylocopa, species representing 22 of 51 subgenera were used (Hurd & Moure, 1963), or 22 of 33 subgenera if Minckley (1998) is followed.These species represent all the major taxonomic divisions in the genus, as recognized from phylogenetic analyses by Minckley (1998), as well as representing all continents on which Xylocopa species are known.Multiple species representatives for the subgenera Koptortosoma and Neoxylocopa were used, because these two subgenera comprise the majority of the described species in Xylocopa .Previous molecular analyses (Leys, Cooper & Schwarz, 2000) has shown that the subgenus Koptortosoma is paraphyletic, showing an African and Oriental clade.Species representatives for Koptortosoma and Neoxylocopa were therefore chosen from a wide geographical range.The majority of subgenera not included in the present study consist of one or a small number of species which, according to Minckley (1998), form well supported clades with the subgenera included here.The set of species studied here was also used in earlier phylogenetic analyses based on sequences of CO1 and Cyt-b (Leys et al ., 2000;GenBank accession numbers AY005222-5275).A list of outgroup and ingroup species are presented in Table 1.GenBank accession numbers of the sequences of EF-1 α are AY005276-5302 and of PEPCK are AY005303-5329.Tissue samples and voucher specimens are kept at the South Australian Museum, Adelaide.

M OLECULAR M ETHODS
DNA extractions: DNA was extracted from flight muscles or legs of bees preserved in 95% ethanol using the CTAB-method (Doyle & Doyle, 1990) with minor modifications.Namely, proteinase K (0.25 mg mL − 1 ) was used during extraction instead of mercaptoethanol.
PCR products were purified using glass-milk in a 'BRESA-CLEAN' kit (Geneworks, Australia) or using columns in a UltraClean™ PCR clean-up DNA purification kit under conditions specified by the manufacturer (MoBio Laboratories Inc.).
PCR products were sequenced using the PRISM Ready Reaction Dye-deoxy terminator cycle sequencing kit (ABI-Perkin Elmer) in 20 or 10 µL reaction volumes according to the manufacturers instructions.PCR primers were used as sequencing primers.Reac- tion products were purified by ethanol precipitation (as specified by ABI) and sequenced on ABI 377 (version 2.1.1,3.0 or 3.3) automated DNA sequencers.Usually two individuals per taxon were sequenced on both strands for both genes to assess intraspecific variation and to check for misidentification, contamination and sequencing errors.Amplified PEPCK fragments varied in length because of two introns of variable sizes.Usually the amplified fragments were too long to obtain a completely readable sequence, but aligning sequences of both strands resulted in perfect matches of approximately 300 bp in the middle of the fragment.Sequences were manually edited and aligned using the program SeqEd version1.03(ABI).

PHYLOGENETIC ANALYSES
Phylogenetic analyses of aligned sequence data were carried out using applications in the program PAUP*4.0b2b (Swofford, 1998), unless otherwise stated.The methods for phylogenetic analyses used in this paper are similar to the methods used for analysing the same set of specimens using mitochondrial genes (see Leys et al., 2000).
Base composition bias for the different codon positions of the two genes were calculated and Chi-square tests were performed to check for taxa with deviations in nucleotide composition.Significantly different nucleotide compositions may be the result of sequencing a copy of the gene with a different evolutionary history from the remaining taxa, but could also be the result of different directional mutation pressures.
Partition homogeneity (ILD) tests (Farris et al., 1995) were performed to evaluate whether incongruence between the two nuclear data sets might be present, but not as a test of combinability (Messenger & McQuire, 1998).Recently, Wenzel & Siddall (1999) recommended that it is always worth combining data sets, because the signal from the separate data sets can be additive, while noise is averaged, therefore creating a 'new' data set with relatively more signal.For this reason, we analysed the data separately when the ILD tests showed significant heterogeneity, but also carried out phylogenetic analyses using combined sequence data from EF-1α, PEPCK and the mtDNA COI and Cyt-b genes (Leys et al., 2000).By using the mitochondrial and nuclear genes in a combined analysis the number of informative characters increases, which might positively influence the resolution of the resulting trees.To investigate which nodes in the phylogeny are supported by the multiple data sets and which nodes have conflicting support from the different data partitions, the program TreeRot (version 2; Sorensen, 1999) was used to calculate partitioned Bremer support indices (Bremer, 1988).To calculate partitioned Bremer support indices heuristic searches with 100 random addition replicates were performed for constraint trees lacking the node of interest (Sorensen, 1999).
Maximum parsimony (MP) analyses were performed using heuristic searches with 100 random additions of sequences and TBR branch swapping to search for the most-parsimonious trees from different islands of trees.Bootstrapping with 1000 pseudoreplicates and the heuristic search option was used to examine the robustness of clades in the resulting trees, and tree lengths (TL), consistency indices (CI) and retention indices (RI) were used as measures of fit.MP analyses were carried out for the nuclear data sets separately, but also in a combined ('total evidence') analysis with mitochondrial (Leys et al., 2000), and morphological (Minckley, 1998) data sets.
Phylogenetic analyses were also carried out using maximum likelihood (ML).To test which model of molecular evolution was the most appropriate for our data sets, a series of nested likelihood ratio tests were performed on the EF-1α and PEPCK data sets separately, using the method described in Huelsenbeck & Crandall (1997).Instead of using a Jukes-Cantor tree as the starting tree, as proposed by Huelsenbeck & Crandall (1997), we used the MP tree with the highest fit scores based on the previous MP analyses.Likelihood scores of this tree, under the assumption of different models of molecular evolution, were calculated using PAUP* and the significance of the likelihood ratios were tested against a Chi-square distribution using the program MODELTEST (version 2.0; Posada & Crandall, 1998).
The large size of our data set, and the large number of parameters that needed to be estimated, made it computationally very intensive to estimate parameters using ML searches.We used a method to estimate the model parameters that reduced calculation times by more than a factor of 30 compared to the above mentioned method.All most-parsimonious trees resulting from a previous analysis were used to calculate likelihood scores and estimate the parameters for the evolutionary model that fitted our data best.The parameter values of the tree with the highest likelihood score were used as starting parameters in the first search of a series of ML searches using the PUZ-ZLE option in PAUP*.PUZZLE in PAUP* is based on the quartet puzzling method of Strimmer & von Haeseler (1996), a fast tree search algorithm using ML methods, which also automatically assigns estimations of support to each internal branch.A number of iterative, heuristic searches with previously estimated parameters, using the resulting tree for estimating new parameters, were performed with PUZZLE until parameters, likelihood scores, and tree topology did not change.These 'final' parameters were used for a ML heuristic search.Our earlier phylogenetic study using mtDNA showed that these estimated parameters were very similar to the parameters estimated by ML with complex models of DNA substitution and a heuristic search (Leys et al., 2000).Due to time constraints ML heuristic searches were done with simple addition of sequences, but without bootstrapping.Maximum likelihood analyses were carried out for the nuclear data sets separately using PAUP* and in combination with the mitochondrial data sets (Leys et al., 2000).
As framework for the analyses of historical biogeography of the carpenter bees a phylogenetic tree was used inferred from Bayesian analysis (Larget & Simon, 1998) using the mitochondrial and nuclear data sets.The advantage of using Bayesian analysis over parsimony or likelihood tree reconstruction methods is that it gives reliable information about the correctness of the obtained tree topology, as well as information on the (posterior) probability of the individual nodes (comparable with bootstrap values), given the data and the model of molecular evolution used (Huelsenbeck et al., 2001;Rokas & McVean, 2000).Moreover, it was estimated that using the Bayesian analysis would reduce calculation times by a factor of at least 20 compared to the time needed for ML bootstrapping.The HKY model of molecular evolution (Hasegawa, Kishino & Yano, 1985) was used with site-specific parameters for seven data partitions in a combined analysis: first and second nucleotide positions of mt-genes combined, and for each of the n-genes; third nucleotide positions of mt-genes combined, and for each of the n-genes; and introns of PEPCK.Analyses were run on a SUN computer following the suggested procedure of the MCMC convergence in the BAMBE (Version 2.01) manual (Simon & Larget, 2000).

ESTIMATION OF DIVERGENCE TIMES AND INFERRING HISTORICAL BIOGEOGRAPHY
Divergence times were estimated using the tree topology resulting from the Bayesian analysis.Likelihood ratio tests using this tree and MODELTEST (version 3.06, Posada & Crandall, 1998) were performed to investigate whether a global molecular clock fitted our data.Non-parametric rate smoothing and penalized likelihood methods in the computer program R8S (Sanderson, 2001) were used to recalculate the branch lengths of the tree according to a relaxed molecular clock with local rates at the branches (Sanderson, 1997).Powell's method for optimizing the objective function (Press et al., 1992) was used and stability of solutions was investigated by rerunning the optimization 20 times with random perturbations of initial solutions (Sanderson, 1997).To calculate divergence times for the nodes in the tree, it was necessary to calibrate the divergence time of at least one node in the tree.This was done by applying the program QDATE (Rambaut & Bromham, 1998) and using data from fossil (Oligocene-Miocene) carpenter bees.The linearized tree with divergence time estimates was then used to infer the time frame in which the diversification of the carpenter bees most likely took place.MacClade 3.08 (Maddison & Maddison, 1998) was used to reconstruct ancestral geographical distributions in the phylogeny, and the atlas of Mesozoic and Cenozoic coastlines of Smith et al. (1994) was consulted to infer the historical biogeography of the carpenter bees.

RESULTS
Sequencing of the EF-1α and PEPCK PCR fragments resulted in 400 and 1051 aligned nucleotide sites, respectively.Table 2 shows the frequencies of the nucleotide bases, nucleotide bias, overall Chi-square values of tests for homogeneity of base frequencies across taxa, and percentage and number of parsimony-informative sites, for the different codons.Base composition bias across taxa was found to be significant only for third codon positions in EF-1α.G-C bias in third codon positions of EF-1α was apparent, varying from 66.4% in the outgroup, Brevineura, and 83.5% in Xylocopoides to 96.2% in Mesotrichia, Hoploxylocopa and Alloxylocopa.Between subgenera, PEPCK sequences varied greatly in length because the amplified fragments contained two introns varying between species from 26 to 499 bp in length.
The two nuclear data sets were shown to be significantly heterogeneous by a partition homogeneity (ILD) test (P = 0.01) applied with invariant characters excluded (Cunningham, 1997).This result is possibly due to the relatively low number of informative sites in EF-1α, most of which were at third codon positions and likely to include considerable numbers of homoplasies, as compared to PEPCK (see Table 2).We therefore analysed the EF-1α and PEPCK data sets separately, but also investigated the effect of combining the data sets for phylogenetic analyses.

MAXIMUM PARSIMONY ANALYSES OF EF-1α
Parsimony analysis of EF-1α resulted in a single island with 350 most-parsimonious trees.The strict consensus tree (Fig. 1a) shows an 'Ethiopian' clade (for consistency this name was adopted from Minckley (1998) although this clade also contains Oriental taxa), supported by bootstrap values of 54-67%, which consists of an Oriental group (including Koptortosoma aruana, K. lieftincki, Platynopoda, Hoploxylocopa and Alloxylocopa, except for the African subgenus Afroxylocopa) and an African group (including Mesotrichia, Koptortosoma pubescens and K. scioensis), and the two subgenera Nyctomelitta and Proxylocopa united with bootstrap values of 74%.Other apparent groups, like the American group, were not supported by bootstrap values larger than 50%.The lack of resolution in the results from the EF-1α data set is most likely due to the low number of parsimony-informative sites (Table 2), with most of the variation in third codon positions.

MAXIMUM LIKELIHOOD ANALYSES OF EF-1α
The same data set as that used for MP analyses was used for ML analyses.The program MODELTEST  1994) with four discrete rate categories is the most suitable model for our data set.Figure 1b shows the ML phylogeny.The topology of the tree is very similar to the trees resulting from the MP analyses.However, the positioning of Lestis within the Ethiopian group and Copoxyla, showing a long branch, breaking up the Neoxylocopa group, has not been found in the previous MP analysis nor in the mtDNA analyses (Leys et al., 2000).

MAXIMUM PARSIMONY ANALYSES OF PEPCK
The PEPCK sequences contained two introns which were variable in length across taxa (subgenera) due to numerous deletions or insertions (indels) of different sizes.Exon/intron junctions were conserved AA(or AG)/GT sequences, and intron/exon junctions AG/AT(or GA) sequences.Aligned intron sequences showed numerous deletions varying in total size from 4.5% (25 sites in Ctenoxylocopa) to 87.1% (499 sites in Brevineura) of the total aligned intron sequence of 573 bp.Indels of similar sizes are shared among (supposedly) closely related taxa and therefore the positions and the different sizes of the deletions potentially can hold phylogenetic information.Alignment gaps in the introns of the PEPCK sequences were coded using the 'simple indel coding' method described in Simmons & Ochoterena (2000).This method is implemented by coding all gaps that have different 5′-and 3′ ends as separate presence/absence characters, regardless of whether such gaps overlap each other or not.Potentially more phylogenetic information may be retained when indels that (partly) overlap are treated as single characters with different character states.However, we could not apply this 'complex indel coding' method as proposed by Simmons & Ochoterena (2000), because the nearly complete lack of intron sequence in some of our taxa (e.g.Brevineura; Biluna) would simultaneously reduce the number of gap characters but also increase the number of its character states.Coding indels according to the 'simple indel coding' method resulted in 54 additional characters of which 18 were parsimonyinformative.
PEPCK data was analysed using the total sequence data (exon and intron sequences, with gaps coded as missing characters) combined with the indel dataset, and using the gap character data set only to explore the phylogenetic information content of the indel characters.Analyses of the combined data set resulted in three resolved clades: the American group (supported by a bootstrap value of 93%) with Nyctomelitta and Proxylocopa as closest relatives; the Xylocopa s.l.group; and the 'Ethiopian' group with all Oriental taxa in a single clade (bootstrap 69%; Fig. 2a).Analysis using the indel data set resulted in a tree topology resolving the American group with Nyctomelitta and Proxylocopa as closest relatives, as well as the Ethiopean group, but the taxa within the Xylocopa s.l. group remained unresolved (18 parsimony-informative characters, tree length = 69, CI = 0.771, RI = 0.830, tree not shown).Considering that this small data set was able to largely resolve the topology, indicates that there was little homoplasy among the indel characters.

MAXIMUM LIKELIHOOD ANALYSES OF PEPCK
Exon and intron data of PEPCK were also used for ML analyses.Gaps were treated as missing characters, because the present ML models do not apply to the evolution of indels.MODELTEST (Posada & Crandall, 1998) performed on the exon data, indicated that the HKY +Γ +Ι model (Hasegawa et al., 1985), correcting for unequal base frequencies, transition/transversion ratio, gamma-distributed rate heterogeneity and proportion of invariable sites, appeared the most suitable model for the exon data set.Analyses were performed for exons only using the HKY +Γ +Ι model, but also for the exon-intron combined data set using the HKY + SSR (site specific rates) model.Using the HKY + SSR model allowed us to partition the data set into different codon positions and introns, to be able to estimate separate site specific rates for every partition.Both ML analyses show a very similar result, which only differs in the placement of Notoxylocopa, Schoenherria and Diaxylocopa, as the sister taxa of the rest of the American clade (Fig. 2b) and some differences in the branching order for the African and Oriental subgenera belonging to the Ethiopian group.The most striking differences between the ML and MP analyses are the findings that the Ethiopian group remains unresolved in the ML tree, while in the parsimony analyses the monophyly of the Ethiopian group was supported by a bootstrap value of 87% (Fig. 2a).In addition, Biluna and Perixylocopa, which were sister taxa of the rest of the Xylocopa subgenera in the MP analyses, are now positioned as a single sister group of Proxylo- copa, Nyctomelitta and the American clade in the ML tree.

ANALYSES WITH COMBINED DATA PARTITIONS
Combined analyses of EF-1α and PEPCK resulted in a topologies very similar to those of the separate PEPCK analyses (Fig. 2a,b).A combined 'total evidence' analysis was carried out by adding sequence data from two mitochondrial genes CO1 and Cyt-b (Leys et al., 2000), and a morphological data set (Minckley, 1998) to the EF-1α and PEPCK nuclear data sets.
MP heuristic searches with 100 random additions of the 'total evidence' data set, with transitions of the third codon positions for the mitochondrial genes excluded, indels in PEPCK introns coded as simple (Simmons & Ochoterena, 2000), and morphological characters equally weighted and coded as unordered, resulted in two islands with four most-parsimonious trees (Fig. 3a).Three major groups, the American group (bootstrap support 95%) forming a clade with Proxylocopa and Nyctomelitta as sister lineages (bootstrap support 77%), the Xylocopa s.l. group and the Ethiopian group (bootstrap support 100%).Biluna was found to be the sister lineage of all other subgenera which form a clade supported by a bootstrap value of 89% (Fig. 4).The phylogenetic relationships among the three major groups and the subgenus Perixylocopa remained unresolved (bootstrap values <50%).Table 3 shows the partitioned Bremer support indices for each node in the phylogeny (Fig. 3a) for the different data partitions in the total evidence analysis.Some conflict between the data partitions is evident for internal nodes of the Xylocopa s.l.(nodes 6 and 7) and Ethiopian clades (node 2).Further, it is apparent that the low resolution among the major clades (nodes 9 and 10) is mainly due to lack of character information rather than to conflicting data.
ML analyses of the mitochondrial and nuclear combined data set is not straightforward, because the different gene partitions fit different models of molecular evolution which also require different parameter settings.The mitochondrial data set showed high A-T bias, probably resulting in saturation of transitions of the third codon positions (Leys et al., 2000), and results from the present paper showed G-C bias in EF- 1α, but no notable nucleotide bias in PEPCK.Based on these observations it would be ideal to analyse the combined data set using different models for each gene partition.This type of analysis is currently unavailable in PAUP*, and other programs (e.g.PAML & Yang, 1999) do not perform tree searches optimally.We therefore carried out a ML analysis using estimated site specific rates for codon positions for every gene partition (including PEPCK introns) in combination with the HKY model using PAUP*.Site specific rates in combination with the general time reversible (GTR, Rodríguez et al., 1990) model was too parameter rich, and consequently too time consuming for our analyses.

Figure 2. (a)
Strict consensus tree of one island of six most-parsimonious trees resulting from MP analyses of PEPCK combined exon and intron sequences and indels coded as separate characters (see text for explanation of gap coding method), using heuristic searches with 100 random addition sequences (tree length = 914, CI = 0.678, RI = 0.691).Bootstrap values (>50%) resulting from analyses with 1000 bootstrap replicates are shown above branches.(b) PEPCK ML phylogram resulting from a heuristic search with ten random addition of sequences and using the HKY +SSR model using exon and intron sequences (-lnL = 5508, ti/tv-ratio = 3.1878, rate parameters codons: introns = 1.371, 1st positions = 0.223, 2nd positions = 0. 1992, 3rd positions = 1.246).3(b) shows the result of a ML analysis using the HKY-SSR model (site specific rates estimated for all codon positions for every gene and for introns of PEPCK, and using an estimate of the transition/ transversion ratio for the whole data set) with transitions of the third codon positions for the mitochondrial genes excluded, and gaps in PEPCK introns coded as missing characters.The topology of the ML tree is very similar to the total evidence MP tree (Fig. 3a).Differences are found only at those branches that were not supported by high bootstrap values in the MP topology: the placement of Perixylocopa, internal topology of the Xylocopa s.l.clade (i.e.placement of Copoxyla) and Ethiopian clade (i.e.placement of Afroxylocopa).Table 4 shows the site-specific rates estimated from the tree topology of Fig. 3(b) using ML for different codon positions and introns of the four gene sequences.

Megaxylocopa
Bayesian analysis was carried out using BAMBE (v2.01;Simon & Larget, 1999).Several initial runs were needed to find the most suitable parameter update values for our data set.Two final runs each with different random seeds, and 2 million cycles of which every tenth tree topology was sampled for analysis, had very similar results.The sampled trees of both runs were therefore summarized together.The most probable tree topology had a posterior probabil-ity of 0.193.The next most probable tree had a posterior probability of 0.089, and only differed with respect of the placement of Perixylocopa as sister taxon of the (American group (Nyctomelitta + Proxylocopa)), as in Figure 3b.The posterior probabilities of the nodes in the most-probable phylogenetic tree inferred from Bayesian analysis are presented in Figure 4.In general, this tree topology only differs from previous total evidence MP (Fig. 3a) and ML analyses (Fig. 3b) with respect to the placement of the major clades and Perix- ylocopa (Fig. 4).All other posterior probabilities for the nodes are high, and therefore this topology should be suitable as the best phylogenetic framework so far for further evolutionary studies.

ESTIMATION OF DIVERGENCE TIMES
We used the tree topology of the Bayesian analysis (Fig. 4) as the preferred tree for biogeographical analysis.
Estimates of divergence times for clades in the phylogeny are needed before inferences on historical distributions and dispersal can be made.A likelihood ratio test using MODELTEST rejected equal rates of molecular evolution at branches in the phylogeny (H 0 : HKY-SSR+clock, H 1 : HKY-SSR, −2ln∆ = 125.1,d.f.= 25, P < 10 −6 ), implying that a global clock rate can not be used to date divergences in the phylogeny.To linearize the branch lengths in the topology of Figure 4 the penalized likelihood method in the program R8S of Sanderson (2001) was used.Evolutionary rates are allowed to vary between adjacent branches within certain limits and are optimized using a smoothing algorithm (see Sanderson, 1998).To calibrate divergence times of the nodes in the linearized tree, the divergence time of at least one node in the tree must have been estimated.A Miocene fossil from China closely resembling extant Nyctomelitta (Zhang, 1990), and fossils from Switzerland (Zeuner & Manning, 1976) most likely related to Xylocopa s.s.(Middle Oligocene) and Copoxyla (Miocene; Hurd & Moure, 1963) were used to estimate the divergence times of nodes (Fig. 4).This was done using QDATE (v1.1;Rambaut & Bromham, 1998), a program that uses quartets of species with fossil dates.The combined sequence data of the four genes were used in the analysis (HKY-model with gamma rate heterogeneity in four rate classes, α∝ = 0.690, ti/tv = 1.054).Using the fossils it was possible to estimate divergence times of two nodes in the phylogeny: the basal node in Xylo- copa s.l.-group estimated using quartets (Xylocopa s.s., Ctenoxylocopa) (Copoxyla, Lestis) revealed a divergence time of 31.1 Mya and node 2 (Fig. 4) using (Proxylocopa, Nyctomelitta) (Xylocopa s.s., Ctenoxylocopa) revealed 39.9 Mya.Table 5 shows the estimated divergence times and 95% confidence limits of the nodes calculated using the program QDATE and the divergence times of other nodes in Figure 4 calibrated with these data using penalized likelihood and rate smoothing methods in R8S.Approximate 95% confidence intervals for the 'R8S' data were calculated using maximum likelihood shape and a cut-off value of 4 (Cutler, 2000).An Eocene-Oligocene (34 Mya) fossil carpenter bee from North America, of which its relationship with contemporaneous American Xylocopa subgenera is unclear (Engel, 2001a), was also used in R8S-analyses.A minimum age of 34 Myr was assigned to each of the nodes leading to the two American lineages in the phylogeny separately: Xylocopoides (node 6) and the American group (node 3), while the other nodes were free to vary.Interestingly, both analyses resulted in similar ages of each of the American lineages.Because the American fossil is slightly older than the other fossils used here, the estimated divergence times in the phylogeny became also slightly older (Table 5).
The inferred ranges of divergence times associated with the nodes 1-4 (Fig. 4, Table 5) show that the evolution of the main subgeneric groups of Xylocopa took place around 55-35 Mya which is well after the initial break-up of Gondwanaland, approximately 148 Mya, and also well after the final break-up of South America and Africa approximately 100 Mya (Smith et al., 1994).

DISCUSSION
Phylogenetic analyses of the carpenter bee subgenera using the nuclear genes EF-1α and PEPCK in combination with the mitochondrial DNA genes Cyt-b and CO1 revealed the same three major clades of taxa that were found in previous analyses based on mtDNA data only (Leys et al., 2000) and analyses based on morphological data (Minckley, 1998).The present study, however, unlike these previous studies, demonstrates higher levels of support for the branches leading to those clades: the American clade together with the subgenera Proxylocopa and Nyctomelitta (bootstrap support 77% MP analysis, 0.866 pp (posterior probability) Bayesian analysis); the Xylocopa s.l.clade (bootstrap support <50% MP analysis, but 0.995 pp Bayesian analysis); and the Ethiopian clade (bootstrap support 100% MP analysis, and 1.0 pp Bayesian analysis; Fig. 4).The monophyly of the subgenera Proxylocopa and Nyctomelitta and the American clade is also supported by high support values (Figs 3a, 4).Further, Biluna is placed basal to all other subgenera included in this study which form a monophyletic group supported by high support values (bootstrap support 89% MP analysis, 0.858 pp Bayesian analysis).In several analyses in the present study based on single gene data, Perixylocopa is found in such a basal position instead (Figs 2a,b and 3a) or shares a clade with Biluna (Fig. 3b).Biluna and Perixylocopa were always found in the same group analyses (Rhysoxylocopa group in Minckley, 1998).Because in our analyses the relationships of the three major clades and Perixylocopa are poorly resolved, Perixylocopa might well share the basal position in the phylogeny with Biluna.However, Bayesian analysis showed that the posterior probability of a tree topology where Biluna and Perixylocopa share a basal clade is less than 0.001, indicating that it is unlikely that both clades are sister lineages.
Overall, we propose that the phylogenetic tree resulting from the combined Bayesian analyses of the four genes EF-1α, PEPCK, Cyt-b and CO1 (Fig. 4), provide the most robust phylogenetic hypothesis so far for large carpenter bees and that therefore the phylogeny should provide a reliable framework for future evolutionary studies.Moreover, the fact that the phylogeny is based on three independent nucleotide data sets (mitochondrial, and nuclear genes), indicates that the tree topology reflects an organismal rather than a gene tree.

DIVERGENCE TIMES
The estimated divergence times inferred from the analyses allowing different rates of molecular evolution give a fair idea of the time in which the major splits in the evolution of the carpenter bees might have occurred.However, it should be noted that the divergences in the linearized tree (Fig. 4) were calibrated using fossil data and although the oldest fossils known for the lineages were used in the calibration, they still represent the minimum ages for these lineages.The actual calculated divergence times are therefore also minimum estimates.It would be difficult to assess the maximum time of divergence in the carpenter bees.One way to get a maximum estimate of a divergence time is to use molecular clocks calibrated using lineages that dispersed to volcanic islands of which the emergence dates are known.Several such data sets exist for the Hawaiian and Canary Islands and several groups of organisms, e.g.drosophilids (DeSalle et al., 1987), beetles (Juan, Oromi & Hewitt, 1995, 1996;Emerson, Oromi & Hewitt, 2000), butterflies (Brunton & Hurst, 1998), but are not available for carpenter bees.However, some idea of a maximum time of divergence comes from consideration of higher order relationships of bees and known fossils.Bees are thought to have evolved from a spheciform wasp ancestor sometimes during the Early Cretaceous (Michener, 1979;2000;Engel, 2001b).The radiation of the major bee groups should have taken place from the Upper to Middle Cretaceous (125-90 Mya;Engel, 2001b), more or less contemporaneously with the major radiation of the flowering plants (Crane, Friis & Pedersen, 1995).By using the higher level phylogeny of the bees (Roig-Alsina & Michener, 1993), in combination with extensive data on amber fossils, Engel (2001b) hypothesized that the family Apidae, of which the genus Xylocopa belongs, appeared around 90 Mya.The subfamily Xylocopinae (including Xylocopa and related tribes) is considered a basal clade in the Apidae (Roig-Alsina & Michener, 1993).Moreover, the oldest single bee fossil, a Trigona (subfamily Apinae), is from Late Cretaeceous (Maastrichtian, 80 Mya) amber of New Jersey (Michener & Grimaldi, 1988).The Apinae possibly diverged later than the Xylocopinae (Roig-Alsina & Michener, 1993), it therefore might be justified to consider that the Xylocopinae arose sometime before 80 Mya.However, these views based on a restricted number of fossils will remain somewhat speculative until a higher level molecular phylogeny of the Apidae, combined with the existing fossil data and application of molecular clock methods, can shed light on this matter.The molecular data presented in this paper clearly shows that the relationships of the extant Xylocopa lineages are narrowed down to 55-35 Mya.

HISTORICAL BIOGEOGRAPHY
The combination of the molecular phylogenetic data, including the estimation of divergence times presented here and the present distributions of the subgenera allow discussion of the origin of the genus Xylocopa.Michener (1979Michener ( , 2000) ) repeatedly discussed the likelihood of alternative hypotheses for the origin (in time and place) of the bees in general, and for families, tribes and genera.The alternative hypotheses are: (1) origin(s) that predated the break-up of Gondwanaland, suggested mainly by disjunct distributions of related bee taxa (e.g.related Meliponini in South-East Asia and South America, or related groups of Paracolletini in South America, Africa and Australia); and (2) origin(s) that postdate Gondwana break-up times.The younger age of several bee groups is suggested by their existence in neighbouring continents mediated by intercontinental dispersal.
Our data, based on the divergence time estimates (<55 Mya) support the hypothesis of a post-Gondwanan origin of the genus.The major splits in the carpenter bee phylogeny appear to have occurred well after the major break-up of Gondwanaland, which ended with the separation of South America and Africa, approximately 100 Mya.The major divergences in the Xylocopa phylogeny most likely took place during the multiple fissions of Laurasia and North America, Africa and the Americas, but before Miocene geographical fusion events of Africa-Europe, India-Asia and North and South America (Smith et al., 1994;Hedges et al., 1996 fig. 6).
The climatic conditions during that period on the major parts of the continents were suitable for Xylo- copa during most of its evolutionary history.Evidence that relatively high latitudes in Europe and Northern America had a tropical climate is provided by Trigona species from Late Eocene Baltic amber (Zeuner & Manning, 1976) and Late Cretaceous amber of New Jersey (Michener & Grimaldi, 1988).Extant Trigona has similar climatic preferences to Xylocopa.
An important factor in the following discussion on the historical biogeography of the carpenter bees is the notion that since the separation of South America and Africa, approximately 100 Mya, all southern continents were in rather isolated positions for nearly 80 Myr.Theoretically the origin of Xylocopa should have been on one of the four continents: South America, North America, Eurasia or Africa.Michener (1979), argued that speculation as to areas of origin was premature using the data available then.However, Engel (2001a) hypothesized, based on higher level bee phylogenies (Roig-Alsina & Michener, 1993), that the most likely origin of Xylocopa would have been Africa.To assess which continent is the most likely origin given our data, ancestral character state reconstructions with MacClade 3.08 (Maddison & Maddison, 1998) and with DIVA 1.1 (Ronquist, 1996(Ronquist, , 1997) ) were performed using distribution areas for the subgenera as character states.Both methods showed that the most likely centre of origin of Xylocopa is Eurasian (Palaearctic or Oriental).This reconstruction is robust with respect to changes of the position of Perixylocopa or to differences in the relative positions of the three main clades, probably because of the fact that the three main groups all contain taxa from Palaearctic or Oriental regions (for distribution data see Table 1).An alternative explanation based on the maximum species diversity of Xylocopa in Africa (Hurd & Moure, 1963), and Michener's (1979) proposed origin of bees or Engel's (2001a) proposal, would be that the origin of Xylocopa is African or South American.This explanation, however, is not supported by the results presented in this paper, and also not by the work of Minckley (1998), both showing an absence of close relationships between African and American subgenera.A possible reason for not finding an American-African connection could be because of insufficient sampling of the African taxa.This explanation, however, seems unlikely because close relationships of the African taxa with Proxylocopa and/or Nyctomelitta, both sister taxa of the American group, were not found in any of Minckley's (1998) analyses, which were based on a nearly complete sampling of subgenera.Sister group relationships of South American with African taxa were also rejected for 14 neotropical long-tongued bee genera (Roig-Alsina & Michener, 1993).Moreover, an African origin for the genus Xylocopa, in combination with our divergence time estimates would imply intercontinental dispersal events over long distances at times (55-35 Mya) when the Southern continents were still well separated, and/or much more recent (<25 Mya) simultaneous dispersal events out of Africa of diverse subgenera from different clades without leaving traces of related subgenera.Both possibilities are considered highly unlikely and even more so in the event of the discovery of older fossil carpenter bees that increase the age of the divergences.
When the carpenter bee phylogeny, estimates of divergence times and the inferred Eurasian centre of origin of Xylocopa (Fig. 5a) are combined with palaeogeographical data (Smith et al., 1994), it becomes possible to reconstruct the historical biogeography of the carpenter bees by assuming a restricted number of independent intercontinental dispersal events (Fig. 5b-f).Although carpenter bees are good fliers, they are considered poor dispersers across large stretches of water (Hurd & Moure, 1963;Hurd, 1978;Michener, 1979).Intercontinental dispersals are therefore expected only to have occurred across land bridges and island chains (e.g.Bering Strait, Caribbean Islands chains), or when continents had drifted sufficiently close together.Based on a Eurasian origin of the genus, the following hypothecial scenario of dispersal events are proposed: (i) strong evidence for a close relationship of East-Palaearctic/Oriental and American taxa is the well supported sister group relationship of Proxylocopa and Nyctomelitta to the American clade (Fig. 4), which suggests that the American lineages evolved from a common ancestor that dispersed into America from Eurasia, most likely across the Bering Strait, sometime before 34 Mya (Table 5, Fig. 5b).This possibility was previously discussed by Hurd & Moure (1963: 22) for the Xylocopa subgenera Mimoxylocopa (East Palaearctic) and Stenoxylocopa (American).Michener (1979) suggested Bering Strait dispersals for several other bee groups.A continuous connection between Asia and North America existed from mid-Cretaceous until the Pliocene (Hamilton, 1983;McKenna, 1983) and has served as a dispersal route for several organisms such as cynipoid wasps (Nordlander, Liu & Ronquist, 1996), mammals and plants (Cox & Moore, 1973: 184-186).(ii) The Xylocopa s.l.group, which comprises several West Palaearctic subgenera, some African genera and the American subgenus Xylocopoides, might have diverged from ancestral Palaearctic Xylocopa by vicariance during fission of the European plates from the rest of Eurasian continent before 37 Mya.The ancestor of the North American Xylocopoides probably dispersed across land bridges at the North Atlantic sometime before 34 Mya (Fig. 5b), although dispersal across the Bering Strait can not ruled out either.(iii) Ancestor(s) of Gnathoxylocopa and Ctenoxylocopa (both from the Xylocopa s.l.group) are likely to have dispersed from Western Europe into Africa when the African continent had drifted sufficiently close to Europe to allow intercontinental dispersals, probably not earlier than 25 Mya (Fig. 5c; Cox & Moore, 1973).(iv) Around the same time the North American ancestors of the American clade were likely to have dispersed into South America (Fig. 5d).(v) Upon the closing of the Middle Eastern gap between Africa and the Asian continent since c. 25 Mya, ancestors of the 'African branch' (Mesotrichia, Afroxylocopa and Koptortosoma) of the Ethiopian group could have dispersed into Africa from South-east Asia (Fig. 5d), either as independent dispersals, or as a single dispersal into Africa, followed by a dispersal back to Eurasia by the Oriental clade around 15 Mya.Finally, two independent dispersals from South-east Asia into Australia would have been facilitated by Australia's approach to South-east Asia around 10 Mya, and by low sea levels during Pleistocene glaciations (Fig. 5e): (vi) Lestis from the Xylocopa s.l. group and (vii) a single dispersal event of a Koptortosoma ancestor from the Oriental branch of the Ethiopian group.The Papua-Australian Koptortosoma are probably monophyletic, because they share at least one unique morphological character, the incomplete intercubitus (Leys, 2000).
The numerous independent dispersals of xylocopine taxa with Northern hemisphere ancestors into southern continents might have been partly driven by the worldwide cooling that started at the end of the Oligocene (Wolfe, 1978; 'latitudinal contraction of the tropics', Hurd & Moure, 1963: 21-22).During the Pleistocene glaciations large parts of the Eurasian continent, with the exception of South-East Asia and the Indonesian Archipelago, would have been unsuitable for carpenter bees (Fig. 5e).The same situation would have occurred in Northern America.It is therefore likely that global cooling at the end of the Tertiary may have influenced dispersal from the north into Africa and South America.The fact that the northern continents were largely unsuitable for carpenter bees might have caused extinction of the sister taxa of those lineages that dispersed to other continents (e.g.Lestis), unless the extant sister taxa were not included in the sampling.Recolonization in northerly directions by temperate climate adapted species (e.g.Xylocopa s.s., Copoxyla, Xylocopoides) occurred in the Holocene (Fig. 5f).
The subgenus Perixylocopa has been left out of the discussion until now because its position in the phy- logeny is the least certain.Minckley's (1998) analyses showed that Perixylocopa belongs to a number of subgenera which he synonymised with Xylomelissa.Unfortunately, a historical explanation is not possible at this stage, because other southern African taxa belonging to that group could not be obtained for DNA analyses.Minckley (1998), however, showed that Xylomelissa belongs to a group that he calls the Rhysoxylocopa clade in which Rhysoxylocopa, Nodula and Maaiana are sister groups of Xylomelissa.It is very likely that the Rhysoxylocopa clade will resolve in the phylogeny as a fourth clade, because the included subgenera appear as a monophyletic group in most of Minckley's (1998) morphological analyses.Rhysoxylo- copa, Nodula and Maaiana have Palaearctic and Oriental affinities and it is therefore likely that an ancestor of Xylomelissa dispersed into Africa.Most other subgenera of carpenter bees that were not included in the analyses conducted here can be assigned to the three main groups that are recognized above.Hence, their inclusion would probably not change the hypothesized pattern of dispersals.There is however, one pair of subgenera (Prosopoxylocopa and Zonohirsuta) considered to be closely related (Hurd & Moure, 1963;Minckley, 1998), that might account for another independent dispersal into Africa.Prosopoxy- locopa is endemic in Madagascar and does not seem to be closely related to any other subgenus from the African continent, but it seems related to the Oriental subgenus Zonohirsuta (Hurd & Moure, 1963;Minckley, 1998).Both subgenera are expected to belong to the Oriental branch of the Ethiopian group.
In conclusion, based on phylogenetic, biogeographical, fossil and palaeo-geographical information it is argued here that the genus Xylocopa most likely had an Oriental-Palaearctic origin and that the present world distribution of the Xylocopa subgenera are mainly the result of dispersal events.The best supporting evidence for this hypothesis is found in the sister group relationship of the Oriental-Palaearctic subgenera Nyctomelitta and Proxylocopa with the Neotropical subgenera.Similar disjunct distributions in other bees have been explained by pre-Gondwanan origins and vicariance.A relevant example is the disjunct distribution of the two meliponine bee (sub)genera Plebeia and Tetragona that both occur in the Oriental or Papua-Australian region and in the Neotropics (Michener, 1979;Roubik, 1989).The similarity with the Oriental-Neotropical affinities in Xylocopa are striking, and therefore a dispersal explanation of the disjunct distributions of meliponines should be considered a strong possibility, especially when considering that meliponines only relatively recently disappeared from most parts of the northern continents because of latitudinal contraction of the tropics in the Neogene.

Figure
Figure 3(b)  shows the result of a ML analysis using the HKY-SSR model (site specific rates estimated for all codon positions for every gene and for introns of PEPCK, and using an estimate of the transition/ transversion ratio for the whole data set) with transitions of the third codon positions for the mitochondrial genes excluded, and gaps in PEPCK introns coded as missing characters.The topology of the ML tree is very similar to the total evidence MP tree (Fig.3a).Differences are found only at those branches that were not supported by high bootstrap values in the MP topology: the placement of Perixylocopa, internal topology of the Xylocopa s.l.clade (i.e.placement of Copoxyla) and Ethiopian clade (i.e.placement of Afroxylocopa).Table4shows the site-specific rates estimated from the tree topology of Fig.3(b) using ML for different codon positions and introns of the four gene sequences.Bayesian analysis was carried out using BAMBE (v2.01;Simon & Larget, 1999).Several initial runs were needed to find the most suitable parameter update values for our data set.Two final runs each with different random seeds, and 2 million cycles of which every tenth tree topology was sampled for analysis, had very similar results.The sampled trees of both runs were therefore summarized together.The most probable tree topology had a posterior probabil-

Figure 5 .
Figure 5. (a-f) Hypothetical scenario of the historical biogeography of Xylocopa plotted on palaeogeographical maps.The distribution of Xylocopa is indicated with black, the direction of dispersal with arrows.See text for further explanation.

Table 1 .
Subgenera and species of carpenter bees examined

Table 2 .
Frequency of nucleotide bases and G-C bias for different codon positions for EF-1a and PEPCK sequences, the probability (P) of the Chi-square test of nucleotide bias among taxa (d.f.= 78), and the percentage and number (N) of parsimony informative characters *gaps coded as 'missing'.Figure 1.(a)Strict consensus tree of one island of 350 most-parsimonious trees resulting from maximum parsimony analyses of EF-1α, using heuristic search with 100 random addition of sequences, and characters unweighted (tree length = 187, CI = 0.597, RI = 0.727).Bootstrap values (>50%) resulting from analyses with 1000 bootstrap replicates are plotted.

Table 3 .
Partitioned Bremer support values for the nodes in Figure 3(a).Values with apparent conflict are in bold

Table 4 .
Site specific rates estimated using maximum likelihood for the different codon positions and introns of four gene sequences (all characters included) Linearized Bayesian analysis tree.Posterior probabilities as values of branch support are shown above the branches.Approximate 95% confidence intervals of the divergence times calculated with R8S are shown with shaded bars.Fossil calibration points are indicated with black dots.Node numbers refer to those in Table5.

Table 5 .
Divergence time estimates and 95% confidence intervals (in brackets) of the main splits in the carpenter bee phylogeny of Figure4.A: ages estimated using fossils from Eurasia.B: ages estimated using a North American fossil fixed at two different nodes to calibrate the other nodes