ore JHR 49: 25-41 (2016) JOURNAL OF | *0eeriewed opevaccets ural AA a torte (>) Hymenoptera http://jhr.pensoft.net The insertional Society of ymenoptersts. RESEARCH Biogeography and demography of an Australian native bee Ceratina australensis (Hymenoptera, Apidae) since the last glacial maximum Rebecca M. Dew!, Sandra M. Rehan?, Michael P. Schwarz! | Biological Sciences, Flinders University of South Australia, GPO Box 2100, Adelaide 5001 2. 46 College Road, Department of Biological Sciences, University of New Hampshire, Durham, NH, USA 03824 Corresponding author: Rebecca M. Dew (dew0009 @flinders.edu.au) Academic editor: /. Neff | Received 9 February 2016 | Accepted 30 March 2016 | Published 28 April 2016 http://zoobank. org/C3 C89FA2-7 DB9-4430-8 924-800 4381 CD84A Citation: Dew RM, Rehan SM, Schwarz MP (2016) Biogeography and demography of an Australian native bee Ceratina australensis (Hymenoptera, Apidae) since the last glacial maximum. Journal of Hymenoptera Research 49: 25-41. doi: 10.3897/JHR.49.8066 Abstract The small carpenter bees, genus Ceratina, are highly diverse, globally distributed, and comprise the sole genus in the tribe Ceratinini. Despite the diversity of the subgenus Neoceratina in the Oriental and Indo- Malayan region, Ceratina (Neoceratina) australensis is the only ceratinine species in Australia. We examine the biogeography and demography of C. australensis using haplotype variation at 677 bp of the barcod- ing region of COI for specimens sampled from four populations within Australia, across Queensland, New South Wales, Victoria and South Australia. There is geographic population structure in haplotypes, suggesting an origin in the northeastern populations, spreading to southern Australia. Bayesian Skyline Plot analyses indicate that population size began to increase approximately 18,000 years ago, roughly corresponding to the end of the last glacial maximum. Population expansion then began to plateau ap- proximately 6,000 years ago, which may correspond to a slowing or plateauing in global temperatures for the current interglacial period. The distribution of C. australensis covers a surprisingly wide range of habitats, ranging from wet subtropical forests though semi-arid scrub to southern temperate coastal dunes. The ability of small carpenter bees to occupy diverse habitats in ever changing climates makes them a key species for understanding native bee diversity and response to climate change. Keywords Climate change, dispersal, DNA barcoding, bayesian skyline plot, haplotype network, population structure Copyright Rebecca M. Dew et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 26 Rebecca M. Dew et al. / Journal of Hymenoptera Research 49: 25-41 (2016) Introduction Molecular studies have greatly increased our understanding of the antiquity of bees and their historical biogeography, especially with respect to centres of origin and sub- sequent dispersal routes (e.g. Fuller et al. 2005; Hines 2008; Chenoweth and Schwarz 2011; Rehan and Schwarz 2015). Other studies using museum collection data have implicated very recent climate change as a possible factor underlying changes in bee abundances (e.g. Cameron et al. 2011; Burkle et al. 2013), but there are surprisingly few studies that have attempted to infer changes in bee abundance beyond the last 200 years (but see Wilson et al. 2014). In the face of likely future climate change, it is important to understand how bees have responded to past climates so that we may better predict future trends. Two recent studies (Groom et al. 2014a; Lépez-Uribe et al. 2014) have used phylogeographic and coalescent Bayesian Skyline Plot analyses to examine changes in bee abundances for tropical halictine (Halictidae) and euglossine (Apidae) bees respectively. Both studies found a strong response to Pleistocene climates, suggesting that these two faunal groups have been impacted by glacial cycles despite their tropical distribution. Small carpenter bees, Ceratina (Apidae: Ceratinini), of the subgenus Zadontomerus in eastern North America also showed a rapid population expansion approximately 20 kya, linked to post-glacial cycles (Shell and Rehan 2016). However, no studies have examined possible impacts of historical climates on bee species spanning temperate to xeric zones, beyond those using museum records. The small carpenter bee genus Ceratina has a nearly global distribution, occurring on all continents except Antarctica (Rehan and Schwarz 2015). This includes recent- ly colonized remote islands in the Southwestern Pacific and southern Indian Ocean, probably via accidental human agency (Rehan et al. 2010a; Groom et al. 2014b). Ceratina originated in Africa in the late Paleocene or early Eocene and showed rapid long distance dispersal events allowing it to eventually colonize the New World by the late Eocene (Rehan and Schwarz 2015). Interestingly, patterns of radiation in this tribe suggest that major long distance dispersal events have been rare and tended to oc- cur more frequently in the early history of this tribe rather than later on, despite there being few geographical barriers to later dispersal events (Rehan and Schwarz 2015). Physical impediments to dispersal, such as water barriers, are actually believed to have decreased in the more recent history of this tribe (Rehan and Schwarz 2015). Ceratina (Neoceratina) is the sister clade to all other Ceratina subgenera and origi- nated from a dispersal from Africa to the Oriental region in the early Eocene, with some species extending into the Palearctic (Rehan and Schwarz 2015). Surprisingly, there is only a single ceratinine species in Australia, Ceratina (Neoceratina) australensis (Perkins, 1912). This species forms the sister lineage to all other C. (Neoceratina) spe- cies, and its stem age is dated to the middle Eocene. Ceratina australensis has become a model species for understanding simple forms of sociality where both solitary and social forms remain in sympatry (e.g. Rehan et al. 2010b, 2011, 2014). Solitary nests comprise about eighty-five percent of the population and are founded by females that Biogeography of an Australian native bee af disperse from their natal nests (Rehan et al. 2014). Dispersal of females could facilitate gene flow between populations across Australia. Michener (1962) recorded Ceratina australensis from subtropical and temperate regions of eastern Queensland and New South Wales but there have been no further studies of its distribution. Here we use haplotype variation at 677 bp of the mitochon- drial ‘barcoding’ region of cytochrome c oxidase subunit 1 (CO1) from 102 specimens of C. australensis, obtained from southern Queensland, mid-New South Wales, north- ern Victoria and southern South Australia, to examine historical demography and geo- graphical patterns in population genetic structure. Based on the dispersal of females from natal nests we predicted that there would be gene flow between populations, influencing the genetic structure and historical demography of the species. Methods Collecting localities for genetic samples Specimens of Ceratina australensis were collected from four populations: (i) Mildura in northwestern Victoria; (ii) West Beach in metropolitan Adelaide, South Australia; (iii) the Cowra region in central New South Wales; and (iv) the Warwick region in south- eastern Queensland (Figure 1). Our four study sites represent a substantial proportion of the geographic and climatic conditions for the species, covering warm temperate forest, semi-arid riverine woodland, mediterranean coastal dunes and central table- lands. GPS coordinates, collection dates and the number of barcoded specimens from each population are shown in Table 1. Nests, predominantly found in dead stems of plants from the genera Cakile (Brassicaceae), Senecio (Asteraceae), Ferula (Apiaceae) and the species Verbena bonariensis L. (Verbenaceae) were collected during early morn- ings or late evenings. This ensured that the bees would be present in the nest and not out foraging. Nests were collected, as this species is rarely observed on flowers (Michener 1962). Nests were stored on ice or at 4 °C until processing. One randomly chosen adult from each nest was stored in 100% ethanol for DNA sequencing. DNA extraction and sequencing techniques One leg of each specimen was incubated overnight in arthropod lysis solution with pro- teinase K. Extractions proceeded using a glass fibre plate and a vacuum manifold to pull the eluates through the membrane, following the procedures detailed in Ivanova et al. (2006). The DNA extract was stored in 50u1 TLE (10mM TRIS, 0.1mM EDTA pHé8). Forward and reverse primers M070 (5'-TCC AAT GCA CTA ATC TGC CAT ATT A-3') and M080 (5'-TAC AGT TGG AAT AGA CGT TGA TAC-3') were used for PCR amplification of a 700 bp region of CO1. We used 27.5 ul reactions with 0.1 yl immolase as the active enzyme, 1 ul of each M070 and M080, 5 ul of MRT Buffer, 15.4 ul water 28 Rebecca M. Dew et al. / Journal of Hymenoptera Research 49: 25-41 (2016) Table 1. Summary of samples collected from each population of C. australensis. Includes GPS coor- dinates, collection dates, total number of specimens sequenced and the number of unique haplotypes recovered. ; . , ‘ Specimens Population Latitude (S) / Longitude (E) Collection Dates Barcoded Haplotypes Cowra, New South Wales 33°52.78' | 148°45.73' October 2015 Ip 8 . ie . 4 June 2013, October 2013, Mildura, Victoria 34°09.25' / 142°09.58 January 2014, April 2014 42 9 Warwick, Queensland 28°12.85' / 152°02.10' January 2010 30 14 West Beach, s : P ' aL COA 34°56.28' / 138°29.95 June 2012, July 2014 19 1 Climate Zones Hot humid summer Warm humid summer, mild winter Hot dry summer, mild winter Hot dry summer, cold winter Warm summer, cold winter Australia Mild/warm summer, cold winter Figure |. Map of Australia with overlaid temperature and humidity climate zones (Bureau of Meteor- ology 2012) showing sampling locations. New South Wales, green; Queensland purple; Victoria, blue; South Australia, yellow. and 5 yl DNA. The PCR cycle began with 10 min of 94 °C. The annealing stage had 5 cycles consisting of 60 s at 94 °C, 90 s at 45 °C and 90 s at 72 °C followed by 35 cycles of 60 s at 94 °C, 90 s at 50 °C and 60 s at 72 °C. Elongation was 10 min at 72 °C with a final 2 min at 25 °C. The PCR products were cleaned using a vacuum plate with 100 yl TLE, with the cleaned products stored in 30 pl TLE. Final forward and reverse sequencing of the cleaned PCR products was performed by the Australian Genome Research Facility. Biogeography of an Australian native bee 29 Alignment and phylogenetic reconstruction Sequences were imported into GENEIOUS v.6.1.6 (Kearse et al. 2012) for editing and alignment. Reverse and forward sequences were combined into a consensus sequence for each sample. As we were comparing individual base pair changes, no ambiguous or unknown base pairs (including at the end of sequences) could be left in the final align- ment. We aimed to maximize both the sequence length and the number of samples. The sequence length was shortened, so that all samples had base pair data covering the same read length without any missing or ambiguous nucleotides, since missing data can lead to spurious results for coalescent analyses (Ho and Shapiro 2011). The final alignment consisted of 102 sequences of 677 bp in length, with 28 unique haplo- types. All edited sequences are submitted to GenBank (accession numbers KR824844- KR824934 and KU664337-KU664347). An undescribed Neoceratina species from the Solomon Islands, Ceratina (Neocer- atina) “Solomons sp.” was included in the alignment as the outgroup to determine the root of the tree. This species has been shown to be phylogenetically distinct to Ceratina australensis (Rehan and Schwarz 2015). The sequences available for this species did not cover the full length of the 677 bp alignment, so the C’ australensis sequences were shortened to 639 bp for this analysis. This sequence attenuation did not remove any of the unique haplotype information present within the C. australensis sequences. ‘The phylogenetic tree was reconstructed in BEAST v.1.8.1 (Drummond and Rambaut 2007) with a Yule Process tree prior. The substitution model HKY+I+G model was identified as the most appropriate based on Akaike information criterion in JMOD- ELTEST (Guindon and Gascuel 2003; Darriba et al. 2012). The analysis ran for 5 x 10’ generations, logged every 1,000 trees, using a fixed mutation rate of 1.0, with all other parameters set to default. The log files were viewed in TRACER v.1.5 to confirm that the posterior had stabilized. A consensus tree was constructed with TREE AN- NOTATER v.1.8.1 with a burn-in of 10,000 trees (i.e. 10 million iterations; Suppl. material 1). The full alignment was then pruned to contain only unique haplotypes (28 se- quences). The analysis was run following the conditions described above. TRACER again confirmed the posterior of the analysis had stabilised and a consensus tree with a burn-in of 10,000 trees was generated (Figure 2). Haplotype networks Analysis of Molecular Variance (AMOVA) was implemented in ARLEQUIN 3.11 (Excoffer et al. 2005) to compare genetic variation within and among Victoria, New South Wales, South Australia and Queensland populations. For these analyses we used all sequences and included all codon positions. All four populations were compared in the full model followed by pair-wise comparisons of each possible pairing in subse- quent AMOVA analyses. The full alignment was imported into NETWORK (Fluxus 30 Rebecca M. Dew et al. / Journal of Hymenoptera Research 49: 25-41 (2016) KK = ME | NTH Se | **K La es | NTH2 KK =] 4 ' ma == eal ae HW) stHi Pa] ry = cd] |_§_ << ——m s! = | Queensland mm — Victoria | = South Australia © — New South Wales mm — =a iy — Figure 2. Maximum credibility tree of the C. australensis haplotypes from Bayesian BEAST analysis. Clades North 1 (NTH1), North 2 (NTH2) and South 1 (STH1) are indicated. Posterior probability values: **--1.0; a 9552-85; Engineering 2016) and a haplotype network was constructed using a median-joining analysis with epsilon set as zero (Bandelt et al. 1999; Figure 3). HAPLOVIEWER confirmed the final network and was used to generate a publication quality figure (Salzburger et al. 2011). Historical demography We used Bayesian Skyline Plot (BSP) analyses implemented in BEAST and TRACER to explore historical changes in effective population size of Ceratina australensis. For these analyses we included all sequences available including duplicate haplotypes, as analyses of only unique haplotypes can give erroneous results (Grant 2015). BSP anal- yses assume that genetic markers evolve neutrally (Ho and Shapiro 2011). The very small number of amino acid changes in our alignment suggests that purifying selection may be operating on 1* and 2™ codon positions, so we restricted BSP analyses to 34 codon positions only. In these analyses we used a GTR model for nucleotide substitu- tions, but did not include an invariant sites parameter (I) since 3 codon positions should not be constrained by selection. Analyses were run for 100 million iterations, sampling every 10,000" iteration to reduce autocorrelation, and were repeated three Biogeography of an Australian native bee 31 NTH2 Figure 3. Haplotype network of C. australensis populations. Each circle represents a unique haplotype, with the numerals inside indicating the number of individuals sampled of that haplotype. Each step be- tween haplotypes indicates one base pair substitution. times to check for convergence. We implemented a strict molecular clock with rate of 1.0, which allows us to readily convert mutations per site per generation into chrono- logical years. We converted the Bayesian Skyline plot scale to chronological years through di- viding it by mutation rate and the number of generations per year. We used the mi- tochondrial mutation rate observed in Drosophila melanogaster Meigen, viz. 6.2 x 10° mutations per site per generation as an estimate of the mutation rate for Ceratina australensis (Haag-Liautard et al. 2008). This method follows a previous study on de- mographic history in Fijian halictine bees in the genus Lasioglossum (Groom et al. 2014a) and North American Ceratina species (Shell and Rehan 2016). We note that the mitochondrial mutation rate for Caenorhabditis elegans (Maupas) is estimated as 9.7 x 10° mutations per site per generation, close to the rate for D. melanogaster, and that they have mitochondrial AT biases of 76% and 82% respectively. Previous stud- ies have reported an AT bias of 74% for the same barcoding region as in our study (Groom et al. 2014a; Shell and Rehan 2016), and our Ceratina haplotypes had an AT bias of 78%. The number of generations was determined as two per year based on nest contents data from the Victorian and South Australian sites (Dew and Rehan, unpub- lished data), which also corresponds to detailed studies on the Queensland population (Rehan et al. 2010b, 2011, 2014). In order to determine whether inferred changes in historical population size in the BSP analysis were significant we also carried out another coalescent analysis using the same parameter settings as our BSP analysis, but implementing a constant population size model. This was then compared to our BSP analysis using a Bayes Factor test. 32 Rebecca M. Dew et al. / Journal of Hymenoptera Research 49: 25-41 (2016) Inferring ancestral distributions Ancestral distributions were inferred using BEAST ancestral traits reconstruction. The full alignment of 102 sequences was used. Sample location for each sequence was coded as a discrete trait (either New South Wales, Queensland, South Australia or Vic- toria). The analysis ran for 2x10° generations, logged every 1,000, with a Yule process tree prior. A HKY+I+G site model with a strict clock of rate 1.0 was employed. All parameters for phylogeny and ancestral trait reconstruction had reached stability, as viewed in TRACER with a burnin of 1x10* generations. Using this burnin a consensus tree was constructed in TREE ANNOTATER. Results Haplotype lineages In total 28 haplotypes of Ceratina australensis were found across the four field sites. The haplotype tree along with posterior probability values for node support from our BEAST analysis is given in Figure 2. This analysis indicates the presence of a clade comprising one New South Wales and two Queensland haplotypes, which is highly supported (PP = 1.0) as sister clade to the remaining haplotypes. Our haplotype net- work analysis (Figure 3) indicates that this clade, which we will refer to as NTH1 (due to their relative northern location), is separated from the common ancestor for the remaining haplotypes by seven nucleotide substitutions, none of which involve amino acid changes. Both the haplotype tree in Figure 2 and the haplotype network in Figure 3 indi- cate geographical structuring of haplotypes. Firstly, the haplotypes in NTH1 were not recovered in any of the 99 specimens genotyped from the more southwestern sampling locations of Victoria and South Australia. Secondly, there was another clade compris- ing five haplotypes that were only recovered from the more northeastern localities of Queensland and New South Wales, and we refer to this clade as NTH2. Lastly, the haplotype, which we refer to as STH 1 (due to its southern location), mostly comprised specimens from South Australia, but also some from Victoria, and it was not repre- sented in any of the Queensland or New South Wales specimens. Two haplotypes are shared between Queensland and New South Wales, with one shared between Queens- land and Victoria. Population structure Pair-wise comparisons among all individuals revealed significantly greater sequence di- vergence between populations than within populations for Victoria to New South Wales and Queensland, and for South Australia to all other populations (Table 2). Queensland Biogeography of an Australian native bee 33 Table 2. Ceratina australensis regional population structure. Diagonal indicates average pairwise differ- ences within populations, number in parentheses indicates total number of sequences for that region; above diagonal are average pairwise differences between populations; below diagonal are pairwise F., values. Significant values (p <0.05) indicated in bold. Population structure Queensland New South Wales Victoria South Australia Queensland 6.88736 (30) 6.5303 (0.49) 6.93968 (<0.0001) 6.5 (<0.0001) New South Wales 0.0598 (0.19) 5.49091 (11) 5.03247 (<0.0001) | 5.09091 (<0.0001) Victoria 0.35579 (<0.0001) | 0.26487 (<0.0001) 2.5331 (42) 3.78571 (<0.0001) South Australia 0.42823 (< 0.0001) | 0.55586 (<0.0001) | 0.58507 (<0.0001) 0 (19) Table 3. Tajima’s D and Fu’s F, tests of neutrality within populations. Segregating sites (S), Tajima’s D score and significance value (D p-value), and Fu's F, value and significance values (F, p-value) are pre- sented. Values in bold are statistically significant (p <0.05). Neutrality tests Victoria South Australia S 19 16 16 0 Tajima’s D 1.36657 0.02304 -1.01646 0 D p-value 0.937 0.558 0.186 1.000 Bus F, -25.15767 -6.57395 -26.633 3.4x10°8 F. p-value 0 0.001 0 1 and New South Wales were not significantly different from one another. These results are mirrored by pairwise F,,. calculations (Table 2), suggesting that all populations ex- cept New South Wales and Queensland are genetically differentiated. There was only one fixed base pair difference identified between any of the populations. This was at base 486, which was a thymine in all Queensland and New South Wales samples but an adenine in all South Australia samples (Victorian samples varied at this base). Tajima’s D value indicated neutral evolution between populations (Table 3). Historical demography Our Bayesian skyline plot (BSP) analysis is summarized in Figure 4 where it is tempo- rally aligned with a graph summarizing two temperature proxies taken from Pahnke et al. (2003). In Figure 4 we have used two x-axis scales, one using mutations per site per generation and the other using years before present, based on a mutation rate of 6.2 x 10° mutations per site and two generations per year. Based on the current best estimate for mutation rate these plots suggest an increase in effective population size beginning approximately 20-18 ka, with a peak rate of increase at about 15-8 kya, and a plateauing after about 6 kya. Our BSP plot shows a long period of stasis from approximately 64-20 kya. This is an artifact of the analysis, where signals prior to the last major effective population size change are lost (Grant et al. 2012; Grant 2015). This period of stasis was trimmed in the final figure to show just the plot from 32 kya. 34 Rebecca M. Dew et al. / Journal of Hymenoptera Research 49: 25-41 (2016) Effective population size 4x103 —s«xt0®—“‘2.3.CO;2 Bandelt H-J, Forster P, Rohl A (1999) Median-joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution 16: 37-48. doi: 10.1093/oxfordjournals. molbev.a026036 Bureau of Meteorology (2012) Australian Climate Averages — Climate Classifications. Govern- ment of Australia. http://www.bom.gov.au/jsp/ncc/climate_averages/climate-classifications/ index.jsp [viewed 29 January 2016] Burkle LA, Marlin JC, Knight TM (2013) Plant-pollinator interactions over 120 years; loss of species, co-occurrence, and function. Science 339: 1611-1615. doi: 10.1126/science.1232728 Byrne M, Yeates DK, Joseph L, Kearney M, Bowler J, Williams AJ, Cooper S, Donnellan SC, Keogh JS, Leys R, Melville J, Murphy DJ, Porch N, Wyrwoll K-H (2008) Birth of a biome: insights into the assembly and maintenance of the Australian arid zone biota. Molecular Ecology 17: 4398-4417. doi: 10.1111/j.1365-294X.2008.03899.x Cameron SA, Lozier JD, Strange JP, Koch JB, Cordes N, Solter LF, Griswold TM (2011) Patterns of widespread decline in North American bumble bees. Proceedings of the National Academy of Sciences of the United States of America 108: 662-667. doi: 10.1073/pnas. 1014743108 Chenoweth LB, Schwarz M (2011) Biogeographical origins and diversification of the exoneu- rine allodapine bees of Australia (Hymenoptera, Apidae). Journal of Biogeography 38: 1471-1483. doi: 10.1111/j.1365-3113.2008.00432.x Darriba D, Taboada GL, Doallo R, Posada D (2012) jModelTest 2: more models, new heuristics and parallel computing. Nature Methods 9: 772. doi: 10.1038/nmeth.2109 Drummond AJ, Rambaut A (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7: 214. doi: 10.1186/1471-2148-7-214 Excofhier L, Laval G, Schneider S (2005) Arlequin ver. 3.0: An integrated software package for population genetics data analysis. Evolutionary Bioinformatics Online 1: 47-50. Biogeography of an Australian native bee 39 Fluxus Engineering (2016) Free Phylogenetic Network Software. http://www. fluxus-engineering. com/sharenet.htm [viewed 29 January 2016] Fuller $, Schwarz M, Tierney S (2005) Phylogenetics of the allodapine bee genus Braunsapis: historical biogeography and long-range dispersal over water. Journal of Biogeography 32: 2135-2144. doi: 10.1111/j.1365-2699.2005.01354.x Grant WS (2015) Problems and cautions with sequence mismatch analysis and Bayesian skyline plots to infer historical demography. Journal of Heredity 106: 333-346. doi: 10.1093/jhered/ esv020 Grant WS, Liu M, Gao T, Yanagimoto T (2012) Limits of Bayesian skyline plot analysis of mtDNA sequences to infer historical demographies in Pacific herring (and other species). Molecular Phylogenetics and Evolution 65: 203-212. doi: 10.1016/j.ympev.2012.06.006 Groom SVC, Ngo HT, Rehan SM, Skelton P, Stevens MI, Schwarz MP (2014b) Multiple recent introductions of apid bees into pacific archipelagos signify potentially large consequences for both agriculture and indigenous ecosystems. Biological Invasions 16: 2293-2302. doi: 10.1007/s10530-014-0664-7 Groom SVC, Stevens MI, Schwarz MP (2014a) Parallel responses of bees to pleistocene climate change in three isolated archipelagos of the southwestern pacific. Proceedings of the Royal Society B 281: 20133293. doi: 10.1098/rspb.2013.3293 Groom SVC, Stevens MI, Schwarz MP (2013) Diversification of Fijian halictine bees: insights into a recent island radiation. Molecular Phylogenetics and Evolution 68: 582-594 doi: 10.1016/j.ympev.2013.04.015 Guindon §, Gascuel O (2003) A simple, fastand accurate method to estimate large phylogenies by maximum likelihood. Systematic Biology 52: 696-704. doi: 10.1080/10635 150390235520 Haag-Liautard C, Coffey N, Houle D, Lynch M, Charlesworth B, Keightley PD (2008) Direct estimation of the mitochondrial DNA mutation rate in Drosophila melanogaster. PLoS Biology 6: e204. doi: 10.137 1/journal.pbio.0060204 Hines H (2008) Historical biogeography, divergence times, and diversification patterns of bumble bees (Hymenoptera: Apidae: Bombus). Systematic Biology 57: 58-75. doi: 10.1080/10635150801898912 Ho SYW, Shapiro B (2011) Skyline-plot methods for estimating demographic history from nucleotide sequences. Molecular Ecology Resources 11: 423-434. doi: 10.1111/j.1755- 0998.2011.02988.x Ivanova NV, Dewaard JR, Hebert PDN (2006) An inexpensive, automation-friendly pro- tocol for recovering high-quality DNA. Molecular Ecology Notes 6: 998-1002. doi: 10.1111/j.1471-8286.2006.01428.x Kearse M, Moir R, Wilson A, Stones-Havas S$, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thierer T, Ashton B, Meintjes P, Drummond A (2012) Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28: 1647-1649. doi: 10.1093/bioinformatics/bts199 Lopez -Uribe MM, Zamudio KR, Cardoso CF, Danforth BN (2014) Climate, physiological tole- rance and sex-biased dispersal shape genetic structure of Neotropical orchid bees. Molecular Ecology 23: 1874-1890. doi: 10.1111/mec.12689 40 Rebecca M. Dew et al. / Journal of Hymenoptera Research 49: 25-41 (2016) Michener CD (1962) The genus Ceratina in Australia with notes on its nests (Hymenoptera: Apoidea). Journal of the Kansas Entomological Society 35: 414-421. Pahnke K, Zahn R, Elderfield H, Schulz M (2003) 340,000-year centennial-scale marine re- cord of Southern Hemisphere climatic oscillation. Science 301: 948-952. doi: 10.1126/ science. 1084451 Perkins RCL (1912) Notes with descriptions of new species, on Aculeate Hymenoptera of the Australian region. Annals and Magazine of Natural History 9: 96-121. doi: 10.1080/00222931208693112 Rehan SM, Chapman TW, Craigie AI, Richards MH, Cooper SJB, Schwarz MP (2010a) Mo- lecular phylogeny of the small carpenter bees (Hymenoptera: Apidae: Ceratinini) indicates early and rapid global dispersal. Molecular Phylogenetics and Evolution 55: 1042-1054. doi: 10.1016/j.ympev.2010.01.011 Rehan SM, Richards MH, Adams M, Schwarz MP (2014) The costs and benefits of sociality in a facultatively social bee. Animal Behaviour 97: 77—85 doi: 10.1016/j.anbehav.2014.08.021 Rehan SM, Richards MH, Schwarz MP (2010b) Social polymorphism in the Australian small carpenter bee, Ceratina (Neoceratina) australensis. Insectes Sociaux 57: 403-412. doi: 10.1007/s00040-010-0097-y Rehan SM, Schwarz MP (2015) A few steps forward and no steps back: long-distance dispersal patterns in small carpenter bees suggest major barriers to back-dispersal. Journal of Bioge- ography 42: 485-494. doi: 10.1111/jbi.12439 Rehan SM, Schwarz MP, Richards MH (2011) Fitness consequences of ecological constraints for the evolution of sociality in an incipiently social bee. Biological Journal of the Linnean Society 103: 57-67. doi: 10.1111/j.1095-8312.2011.01642.x Salzburger W, Ewing GB, von Haesler A (2011) The performance of phylogenetic algorithms in estimating haplotype genealogies with migration. Molecular Ecology 20: 1952-1963. doi: 10.1111/j.1365-294X.2011.05066.x Shell WA, Rehan SM (2016) Recent and rapid diversification of the small carpenter bees in eastern North America. Biological Journal of the Linnean Society 117: 633-645. doi: 10.1111/bij.12692 Wilson JS, Carril OM, Sipes SD (2014) Revisiting the great American biotic interchange through analyses of amphitropical bees. Ecography 37: 791-796. doi: 10.1111/ecog.00663 Biogeography of an Australian native bee 4] Supplementary material | Cladogram including Ceratina (Neoceratina) Solomons sp. as the outgroup to root the tree Authors: Rebecca M. Dew, Sandra M. Rehan, Michael P. Schwarz Data type: EPS file Explanation note: Cladogram including Ceratina (Neoceratina) Solomons sp. as the outgroup to root the tree. Posterior probability values: *** = 1.0; ** = 95; * = 85. Copyright notice: This dataset is made available under the Open Database License (http://opendatacommons.org/licenses/odbl/1.0/). The Open Database License (ODDbL) is a license agreement intended to allow users to freely share, modify, and use this Dataset while maintaining this same freedom for others, provided that the original source and author(s) are credited.