The impacts of glacial cycles on the geographical distribution and size of populations have been explored for numerous terrestrial and marine taxa. However, most studies have focused on high latitudes, with only a few focused on the response of biota to the last glacial maximum (LGM) in equatorial regions. Here, we examine how population sizes of key bee fauna in the southwest Pacific archipelagos of Fiji, Vanuatu and Samoa have fluctuated over the Quaternary. We show that all three island faunas suffered massive population declines, roughly corresponding in time to the LGM, followed by rapid expansion post-LGM. Our data therefore suggest that Pleistocene climate change has had major impacts across a very broad tropical region. While other studies indicate widespread Holarctic effects of the LGM, our data suggest a much wider range of latitudes, extending to the tropics, where these climate change repercussions were important. As key pollinators, the inferred changes in these bee faunas may have been critical in the development of the diverse Pacific island flora. The magnitude of these responses indicates future climate change scenarios may have alarming consequences for Pacific island systems involving pollinator-dependent plant communities and agricultural crops.
Evidence of increases in global average air and ocean temperatures, and rising global sea level, is widely accepted in the scientific community . Such changes are expected to have varying (and potentially very large) effects on ecosystem composition and function . Potential impacts of climate change on pollinators are important because of their role in angiosperm sexual reproduction . Differing phenological responses between plants and their pollinators to changing climates are particularly concerning , yet few empirical studies exist ( and references therein).
Short-term responses for some biota to recent climate change can be assessed using museum/herbaria records extending back over the last few hundred years . However, two problems arise when assessing the impact of climate change on plant–pollinator relationships using data that only extend back over recent history: (i) it can be difficult to discriminate between climate change effects per se and anthropogenic effects unrelated to climate, such as habitat destruction, introduction of exotics and pesticide use ; and (ii) important changes to plant–pollinator relationships may arise over much longer time frames , or short-term disruptions may be reversed as the affected elements adapt to changed conditions.
One way to explore the effects of climate change on plant–pollinator relationships is to reconstruct phylogenies of the various elements, and then to explore diversification patterns and past population sizes. This approach has been used previously to reconstruct how climate fluctuations impact the population dynamics of both terrestrial and marine organisms, including mammals  and fish . Studies of insect populations have also revealed demographic responses to glacial cycles . However, studies of terrestrial impacts of glacial cycles in tropical islands are scarce, despite their role in understanding island biogeography .
Recently, Groom et al.  examined population size and haplotype accumulation of halictine bees in the genus Lasioglossum subgenus Homalictus (herein referred to simply as Homalictus) in the Fijian archipelago. They argued that there was a decline in population size in the Late Pleistocene, approximately 20 ka, followed by a dramatic increase in population size shortly after. They argued that this pattern of decline followed by a rapid bounce-back could be a response to changing climate, with initial contraction of populations during cooling/drying conditions, followed by rapid expansion as global climate reverted to an interglacial period. The effects of human transformation of ecosystems were excluded as a confounding effect in their inferred bee demographics, because the major increase in historical population size occurred before human occupation of the Fijian archipelago (approx. 3000 years ago). However, there might be other possibilities, unrelated to climate, which could explain the observed patterns (e.g. introductions of foreign plant species that provide different food sources, adaptations unrelated to changing climate or pollinator extinctions leading to filling an available niche). One way to further examine whether climate change really did impact on bee population sizes in tropical islands is to examine changes in bee demographics in other isolated archipelagos of the southwest Pacific (SWP). For such a scenario to be supported, changes in populations should coincide with the onset of the last glacial maximum (LGM) at 23 ka and with periods of subsequent warming commencing approximately 19 ka .
The SWP has a very depauperate bee fauna , and most (if not all) of the long-tongued bee species in Fiji represent anthropogenic dispersals [16,17]. The short-tongued Homalictus are the most abundant bee element of many of the SWP islands. In Fiji, there are four Homalictus species [17,18], six species for the archipelago of Vanuatu , and Samoa has the greatest recorded diversity in the SWP with 10 species . Here, we use the mitochondrial cytochrome oxidase c subunit I (COI) gene to examine changes in past population sizes of Homalictus in three isolated (each separated by at least approx. 800 km) island archipelagos; Vanuatu, Fiji and Samoa. We explore the timing of bee colonization in each archipelago and use coalescent analyses to determine whether historical demographic changes show patterns that might indicate common responses to historical climate changes in the SWP.
2. Material and methods
(a) Study system
We collected Homalictus specimens visiting flowers from three isolated neighbouring archipelagos in the SWP: Vanuatu, Fiji and Samoa, spanning a total distance of approximately 4000 km (electronic supplementary material, figure S1). These island systems form highly fragmented groups of complex and relatively recent geological histories, with the largest subaerial land masses for Vanuatu, Fiji and Samoa emerging at approximately 3–8 , 12.5–29  and 5 Ma , respectively. Floral hosts included numerous cosmopolitan and invasive species throughout village gardens, roadsides and native vegetation, with many flowering year-round. Identifications were based on the species descriptions of Perkins & Cheesman , Pauly & Villemant  and Michener .
(b) DNA processing
Extraction, amplification and sequencing of mtDNA was performed by the Canadian Centre for DNA Barcoding at the Biodiversity Institute of Ontario, using standard protocols . Bidirectional sequencing using the universal primer pair LepF1/LepR2  from a single leg of each specimen produced a 654 bp sequence of COI for 425 Homalictus specimens from Fiji, 177 from Vanuatu and 143 from Samoa. All trace files were examined using Geneious Pro v. 6.1.6 , and haplotypes with one or more ambiguous nucleotide in both forward and reverse directions at the same site were removed from analyses. Identical sequences were also removed, leaving 83 unique haplotypes from Fiji, 30 from Vanuatu and 41 from Samoa. In total, 207 unique haplotypes were aligned, with the 154 Pacific representatives combined with 53 outgroup Australian species of the closely related Lasioglossum subgenus Chilalictus and further Homalictus sequences obtained via GenBank (see the electronic supplementary material, table S1). The resulting dataset (base frequencies—A: 33.7%, C: 14.3%, G: 11.4%, T: 40.6%) exhibits strong AT bias (74.3%) at the third codon position, as is common in Hymenoptera . To avoid sequence contamination from potential Wolbachia infection, sequences were screened using the NCBI BLASTn database. All specimens are stored in the Schwarz Bee Collection at Flinders University, South Australia.
Phylogenetic analyses were based on unique COI haplotypes rather than on species, and correspondence between the two is discussed in the Results and Discussion sections.
(c) Phylogenetic analyses
Phylogenetic analyses were performed using Bayesian inference implemented in BEAST v. 1.7.5  with an uncorrelated lognormal relaxed molecular clock for 50 million iterations, sampling every 1000th iteration, under a Yule process speciation model. Chilalictus was constrained as outgroup with a mean root height of 20 Ma (parameter treeModel.rootHeight normal; mean = 20, s.d. = 3.0) as per the multi-gene phylogeny of Gibbs et al. . The most suitable substitution prior was determined as GTR + I + Γ by an AIC criterion in jModelTest v. 2.1.4 , and individual uncorrelated lognormal rate priors set with mean (parameter ucld.mean) following an exponential distribution (mean = 0.33, offset = 0) and standard deviation (parameter ucld.stdev) followed by an exponential distribution (mean = 0.33, offset = 0). Subsequent log files were then examined in Tracer v. 1.5  to determine whether effective sample sizes were adequate (i.e. more than 100). The resulting haplotype tree (electronic supplementary material, figure S2) was compared with an additional genetic-distance-based analysis (electronic supplementary material, figure S3). Our genetic distance analysis used a neighbour-joining technique applied to uncorrected p-distances in PAUP* v. 4.0b . Missing gene fragments were not included when calculating pairwise distances, and trees were constructed using a heuristic search with 100 random taxon–addition sequence replicates, tree–bisection–reconnection branch swapping, and no limit on the number of shortest trees retained.
(d) Patterns of haplotype diversification
We explored temporal patterns of haplotype accumulation using log lineage through time (LTT) plots based only on unique haplotypes for each of the three major archipelago clades. One clade representing a second dispersal into Samoa and two other clades representing later dispersals into Vanuatu (electronic supplementary material, figure S2) were not included in the analyses because of the small number of haplotypes recovered from them. Separate datasets representing each archipelago were then assessed for appropriate models using jModelTest v. 2.1.4 , and assigned a GTR + I + Γ substitution model. Datasets were unpartitioned for codon position, with default priors for gamma distributions, and assigned a constant rate strict molecular clock of 1.0. Fixing the rate to 1.0, rather than using a relaxed clock, allows the resulting branch lengths to be interpreted in units of mean substitutions per site, which can be converted to years given a mutation rate per generation and the number of generations per year. Based on a Bayesian skyline coalescent model implemented in BEAST v. 1.7.5 , analyses were then run for 50 million generations with parameters logged every 5000 iterations. We then produced LTT plots separately in Tracer v. 1.5  for each archipelago.
(e) Estimating changes in Ne
To determine whether the observed changes in haplotype accumulation corresponded with shifts in the effective population sizes (Ne) for the archipelagos, we produced a Bayesian skyline plot (BSP) for each clade based on the same trees produced for our LTT analyses using a constant-rate fixed strict molecular clock of 1.0 in Tracer v. 1.5 . This produces a tree where branch lengths are proportional to the number of substitutions, and is required to be able to determine node ages based on an estimated number of generations per year. Because our dataset includes representatives from across a wide elevation range (0–1100 m.a.s.l.), it is possible that populations at higher elevations have fewer generations per year, and it should be considered that age estimates for these lineages may be underestimated compared with lower elevation populations.
(f) Dating estimates
There are no documented bee fossils in the Oceania region, and our analyses covered recent divergences; thus, calibrating chronogram age estimates is difficult. However, there are other approaches to molecular dating. Few data exist to assess variation in de novo mutation rates among taxa . However, two studies have estimated single-nucleotide mutation rates for invertebrate mitochondrial genomes, on the nematode Caenorhabditis elegans  and the fruit fly Drosophila melanogaster . These respective rates were within a single order of magnitude of each other (9.7 × 10−8 and 6.2 × 10−8 per site per generation), and both taxa had AT biases similar to our Homalictus dataset (76% and 82%, respectively, compared with 74.3% in Homalictus). Consequently, a mutation rate in the region of 10−7 to 10−8 single-nucleotide mutations per generation might be expected for Homalictus. Other estimated rates of mutation in COI of insects were incorporated into divergence date estimates for the bee genus Bombus, where Duennes et al.  applied a rate range of 1.0 × 10−7 to 1.0 × 10−9 substitutions per site per year (one generation per year in Bombus). While these estimates differ across two orders of magnitude, they are not adjusted for varying base composition, which can have very strong effects . Consequently, we deemed the Drosophila rate of 6.2 × 10−8 per site per generation the most appropriate estimate for our study.
Converting the relative age of nodes based on mutation rates into chronological years requires that we know the number of generations per year for Homalictus in our study sites. Unfortunately, the only studies of Homalictus nesting biology are those of Knerer & Schwarz  on a temperate species from southern Australia which produces two generations per year. Studies on temperate Australian Lasioglossum species in the Lasioglossum subgenus Chilalictus also indicate two generations per year . Although in the tropics, Michener & Bennett  demonstrate almost year-round activity in populations of Halictus ligatus, which produced three broods per year. Furthermore, the neotropical halictine Megalopta genalis has been shown to have an egg-to-adult period of approximately 35 days, which would allow for up to four generations per year . The number of generations per year in the tropical habitats of the SWP is therefore likely to comprise three as a minimum, but could be as many as four or more given short egg-to-adult periods and year-round activity.
Using the chronogram produced under the Bayesian skyline coalescent model analysis above, we were able to estimate (i) the likely time at which Homalictus arrived in each archipelago, and (ii) the changes in haplotype accumulation and change in effective population size based on key nodes of the tree topology. By setting the number of generations per year, we are able to convert the age of the nodes provided in mutations per site per generation to calendar years.
Our Bayesian inference analysis indicated a single colonization of Fiji by Homalictus, two colonizations of Samoa and three for Vanuatu (electronic supplementary material, figure S2). The neighbour-joining analysis (electronic supplementary material, figure S3) shows a highly similar topology with clearly monophyletic clades for each of the three Pacific archipelagos, although their arrangement did not agree with our Bayesian approach. The Samoan clade was recovered as sister to the Fijian plus Vanuatu clades. Importantly, the distribution of branch lengths within and between major clades in each archipelago is very similar for these very different tree-building approaches, indicating that both topological and diversification patterns in our analyses are robust to very different analytical approaches.
(a) Haplotype accumulation
LTT plots are given for each of the three archipelago clades in figure 1a–c and are juxtaposed with the corresponding chronograms from our BEAST analysis (figure 1g; electronic supplementary material, figure S4). Although the number of haplotypes differs within each clade, all three curves share a highly similar shape that increases steadily before rapidly increasing at 5 × 10−3 mutations per site (electronic supplementary material, figure S5). The timing of these rapid increases begins at approximately the same time point, which coincides with the emergence of a hyper-diverse haplotype clade in the Fijian representatives and the most diverse haplotype clades in the Vanuatu and Samoan clades (figure 1g). Considered separately, the changes in haplotype accumulation are striking, and together they suggest parallel temporal responses to a shared stimulus able to influence species over large spatial scales.
(b) Effective population size
Modelling the effective population size over time using BSPs (figure 1d–f; electronic supplementary material, figure S6) reveals almost identical patterns of change. The y-axis of these plots can be interpreted as Ne mutation rate per site, and if mutation rates are assumed to be constant, then this can be interpreted as a measure of relative Ne over time. Again, we see largely stable populations over the majority of time, but a dramatic drop in Ne occurs at approximately 5.0 × 10−3 mutations per site before a large increase at 2.5 × 10−3 mutations per site. Aligning the scales of our BSPs and chronograms indicates these changes appear to occur within a very similar time frame, and the tree topologies (figure 1g) demonstrate that this timing corresponds to both a sudden decrease and then sudden increase in Ne across each of the archipelagos.
(c) Age estimates
We were able to determine the approximate ages of two key components of our plots by applying the single-nucleotide mutation rate , and our estimated number of generations per year. First, based on node ages of 0.1322, 0.0999 and 0.051 mutations per site and four generations per year, we obtain crown age estimates of 533, 403 and 206 ka for Vanuatu, Fiji and Samoa, respectively. The crown age refers to the age of the most recent common ancestor for a particular group of taxa. If we instead assume three generations per year, then the estimates increase to 711, 537, 274 ka for Vanuatu, Fiji and Samoa, respectively. These various ages are only a fraction of the islands respective geological ages of 3–8 , 12.5–29  and 5 Ma . Second, by aligning our resulting tree topologies with both LTT and BSPs, we determined that both the increase in haplotype accumulation and change in Ne occurred at approximately 0.006 mutations per site. Based on four generations per year, this equates to 24 ka, or approximately 32 ka if three generations per year are assumed. This suggests a stimulus or stimuli common across all three isolated island groups. Likewise, if the increase in Ne occurred at 0.0025 mutations per site, then the estimated onset of increasing Ne is approximately 10 ka based on four generations per year, or 13 ka based on three generations. Based on these estimates, the changes presented in Ne appear broadly concordant with a major global change in the last glacial termination (figure 1h, with Heinrich event 1 labelled H1).
Numerous studies have outlined the effects of glacial cycles on Palearctic and Nearctic biota, with contractions to southern refugia during glacial maxima and then subsequent northern expansions following deglaciation . Less detailed studies have inferred similar patterns in the Southern Hemisphere , but few have examined how biota in tropical latitudes have responded to these global cycles [44,45]. The tropical Pacific represents a major driving force behind current global climates, and is expected to have played a considerable role throughout glacial periods . Despite its importance, reconstruction of climatic variation in the region over time has been limited by short or insufficiently resolved records, and confidence in extrapolating data from neighbouring regions is therefore similarly limited.
We found that the principal bee fauna in three major SWP archipelagos, spanning a distance of approximately 4000 km and representing multiple colonization events but with no subsequent interchanges, show remarkably similar patterns of haplotype accumulation since the late Pleistocene. Our LTT plots suggest marked anti-sigmoidal (S-shaped) patterns for haplotype accumulation in the Homalictus fauna from Fiji, Vanuatu and Samoa, and our phylogenetic analyses indicate that these similarities cannot be explained by dispersal events linking the bee faunas in these regions. Consequently, Homalictus bees in these regions appear to have independently responded to some factor(s) influencing each of the archipelagos, and any such factor(s) seem to have operated within a close time frame for each archipelago. For any such factor(s) to have impacted over such a large geographical scale, there are few likely explanations.
Anti-sigmoidal LTT plots can arise from both major extinction events , and three-phase Yule models of radiation where cladogenesis is initially high, followed by a period of low speciation and extinction rates, with a second subsequent period of high diversification rates. Discriminating between these alternatives can be very difficult . However, our BSP analyses clearly indicate major declines in effective population sizes of Homalictus for each of the three archipelagos at approximately the same time point. Although single-locus demographic reconstructions are limited in their confidence with increasing time , very recent time scales ensure strong phylogenetic signal for conspecific variation and render the impact of saturation on our analyses negligible. Instead, confidence in estimating this time point is limited owing to potential variation in the applied mutation rate and generations per year, yet all likely scenarios commonly indicate the decrease and subsequent increase occurred either side of the LGM. The LGM involved both global cooling as well as considerable decreases in precipitation for widely separated regions neighbouring the SWP, including Australia , Borneo  and Chile . Such conditions have been shown to have influenced similar contraction/expansion population dynamics in the goby fish of the SWP owing to sea-level fluctuations, determined via deep-sea sediment cores . However, the impacts on terrestrial biota on this region have not been well explored.
A recent study of Borneo stalagmite δ18O records provides the closest comparison for the possible SWP climate variation and demonstrates a 100 kyr glacial–interglacial cycle throughout the Quaternary with a clear increase in atmospheric moisture since the LGM . Palynology studies of the highland Lake Tagamaucia in Fiji  and Lynch's Crater in northeast Australia  indicate a distinct change in vegetation composition beginning approximately 15 ka during a marked post-LGM warming period. With these exceptions, very few studies have attempted to reconstruct SWP climates across the LGM, so the impact of climate change in this region remains largely unexplored.
Our results suggest that some factor(s) forced near-simultaneous changes in bee populations across a wide region of the SWP. We are unable to say whether these changes are due to climate directly forcing changes in bee populations, indirect effects arising from changes in angiosperm communities, or a mixture of both. However, because bees are key pollinators in terrestrial ecosystems, it seems very likely that our inferred fluctuations in bee demography will have had very broad impacts. Island systems often harbour low species richness, which can result in a broadened range of host plant species for pollinators owing to reduced competition, resulting in keystone ‘super generalist’ species . This seems very likely for Fiji, Vanuatu and Samoa, where most (if not all) non-Homalictus bees are likely to represent post-LGM and/or anthropogenic arrivals [16,17].
If the changes in bee populations that we infer here involve corresponding changes in angiosperm communities, either via shared forcing events or effects on pollinator abundance, then events surrounding the LGM are likely to have broadly impacted tropical island ecosystems as well as the better-understood Palearctic and Nearctic biota. Given predicted climate change scenarios, studies are needed to determine how climate change may impact pollinator-dependent plant communities, including agricultural crops. As tropical ecosystems are major sources of biodiversity, incorporating a large proportion of the developing world, the consequences of disrupting plant pollinator interactions could be considerable.
DNA sequences may be accessed via corresponding GenBank Accession or BOLD ID numbers listed for all included specimens in the electronic supplementary material, table S1.
We thank all members of the South Pacific Regional Herbarium at the University of the South Pacific, in particular Marika Tuiwawa and Alivereti Naikatini, for their invaluable assistance with Fijian field logistics and expertise. Vanuatu collections would not have been possible without the assistance of Linette Berukilukilu and Plant Health and Quarantine, who oversaw remote field collections and facilitated permit acquisition in Vanuatu. The support of Afele Faiilagi and the Ministry of Natural Resources and Environment enabled sampling to be conducted in Samoa. This research was possible thanks to the financial support of the Australia Pacific Science Foundation, Rufford Foundation, National Climate Change Adaptation Research Facility, and an Australia Awards Endeavour Research Fellowship awarded to S.V.C.G.
- Received December 17, 2013.
- Accepted April 8, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.