Thrice Out of Africa: Ancient and Recent Expansions of the Honey Bee, Apis mellifera

See allHide authors and affiliations

Science  27 Oct 2006:
Vol. 314, Issue 5799, pp. 642-645
DOI: 10.1126/science.1132772

This article has a correction. Please see:


We characterized Apis mellifera in both native and introduced ranges using 1136 single-nucleotide polymorphisms genotyped in 341 individuals. Our results indicate that A. mellifera originated in Africa and expanded into Eurasia at least twice, resulting in populations in eastern and western Europe that are geographically close but genetically distant. A third expansion in the New World has involved the near-replacement of previously introduced “European” honey bees by descendants of more recently introduced A. m. scutellata (“African” or “killer” bees). Our analyses of spatial transects and temporal series in the New World revealed differential replacement of alleles derived from eastern versus western Europe, with admixture evident in all individuals.

The long-standing association between humans and honey bees, Apis mellifera, is evidenced by 7000-year-old cave paintings depicting honey collection from wild bee nests (1). Managed honey bee colonies satisfy the pollination requirements of both modern agriculture and the demand for products such as honey, wax, and royal jelly (2). With the publication of the honey bee genome sequence (3) and development of high-density single-nucleotide polymorphism (SNP) markers, tools are now available to study evolutionary processes occurring in both native and introduced populations, including Africanization in the New World, and the genomewide consequences of ancient and recent evolution in this important social insect.

The genus Apis comprises 10 species, 9 of which are confined to Asia (4). A. mellifera, however, is distributed from sub-Saharan Africa to central Asia and northern Europe and is composed of more than two dozen morphologically and geographically distinct subspecies (5). It is commonly believed that A. mellifera split from its closest relative, A. cerana, in western or central Asia (where these species' ranges are closest) and subsequently expanded into Europe and Africa (5, 6). Existing populations in Europe and Africa are hypothesized to derive from populations currently found in Asia, where A. mellifera occurs as far east as Kazakhstan (6). However, an apparent discrepancy between the age of divergence among A. mellifera subspecies [0.7 to 1.3 million years (4, 5, 7)] and the split 6 to 8 million years ago (4, 7, 8) between A. mellifera and A. cerana suggests alternative scenarios. One hypothesis [suggested by E. O. Wilson, citing a personal correspondence from C. D. Michener (9)] is that A. mellifera originated in the tropics or subtropics of Africa.

In North America, introductions of the subspecies Apis mellifera mellifera began as early as 1622. Subsequent introductions of A. m. ligustica (the “Italian” bee) began in 1859, followed by introductions of at least seven other subspecies from Europe, the Near East, and northern Africa (the descendants of which are collectively called “European”) (10). Early introductions in South America are less clear, but probably also involved early introduction of western European subspecies (A. m. mellifera and A. m. iberiensis) and later introduction of eastern European subspecies (10). In 1956, a subspecies from the African savannas, A. m. scutellata, was intentionally introduced to Brazil. A. m. scutellata became established and dispersed northward and southward from Brazil through South and Central America (1113), hybridizing with and displacing previously introduced honey bees. Africanized bees reached their southern limit in Argentina in the 1970s (14), but the eventual northward limit in the United States [first invaded in 1990 (15)] is unknown. Although ample evidence shows that both European and African alleles occur in Africanized populations (14, 1621), A. m. scutellata–derived characteristics, including nesting biology, swarming and absconding behavior, foraging, diet, and mitochondrial DNA (mtDNA), tend to displace European characteristics over time (13).

We analyzed 1136 SNP genotypes (22, 23) from 328 A. mellifera, including 175 individuals from native populations (representing 14 subspecies in Europe, Africa, and Asia) and 153 individuals from introduced populations in North and South America. We inferred the ancestral genotype for a subset of SNPs by analyzing 13 additional individuals from 3 congeneric species, A. cerana, A. dorsata, and A. florea.

In the Old World, there was clear population structure at the level of major groups (each containing multiple subspecies) and geographical subspecies. Principal component analysis (PCA) revealed four major clusters (Fig. 1A) consistent with Ruttner's four morphometrically defined lineages (5, 24): M (western and northern Europe), C (eastern Europe), O (Near East and central Asia), and A (Africa). The four groups were broadly consistent with mitochondrial data (25) and with assignment by the program Structure (Fig. 1C) (22). Pairwise FST between major groups ranged from 0.242 to 0.565 (table S2).

Fig. 1.

Patterns of genetic variation in A. mellifera. (A and B) Principal component plots of Old World individuals based on the 1136 SNP loci (22). (A) Bees cluster in four main groupings (M, C, O, and A) in the coordinates of PC1 and PC2. (B) Higher order PCs partition geographical subspecies. See figs. S1 and S2 for additional PCA results. (C) Results of Structure analysis (22). Individuals are represented by vertical lines, grouped by geographical subspecies (Old World) or sampling regime (New World). Division of individuals into colored segments represents the probability of assignment of that individual to each of K groups. K represents an arbitrary number of hypothetical populations, shown on the right. Colors for K = 4 correspond to colors used in PCA (A and B) and Fig. 2, A and C to I.

Each of the 10 geographical subspecies with 9 to 20 individuals was genetically distinct, and subspecies could be partitioned in either PCA or Structure analysis (Fig. 1B; figs. S2 and S3). FST among subspecies was 0.501; much lower values were obtained for subspecies within each of the major groups (0.053, 0.129, 0.118, and 0.082 for M, A, O, and C, respectively). Despite considerable differentiation, admixture between geographically proximal populations in the Old World was evident in both PCA (Fig. 1A) and Structure analysis (Figs. 1C and 2A). Comparison of linkage disequilibrium (LD) between subspecies (fig. S4) also suggested differences related to admixture, bottlenecks, or both. For example, A. m. mellifera and A. m. ligustica exhibited relatively high LD for SNP pairs separated by ≤10 kb (mean correlation coefficient r2 of 0.45 and 0.30, respectively, decreasing to 0.18 and 0.06 for SNP pairs separated by 50 to 100 kb). In contrast, A. m. scutellata and A. m. intermissa exhibited relatively little or no LD (0.08 and 0.13 for ≤ 10 kb; 0.07 and 0.10 for 50 to 100 kb). LD was low overall and affected only a small number of SNPs (fig. S4), consistent with extraordinarily high recombination rates in A. mellifera [∼19 cM/Mb (26)].

Fig. 2.

Geographical and temporal patterns of diversification. (A) Assignment probabilities from Structure analysis (Fig. 1C, K = 4) of Old World individuals. A. m. pomonella (N = 3) collected from Kyrgyzstan. Three A. m. carnica (indicated by asterisk) were provided by a research facility. (B) Neighbor-joining tree based on allele-sharing distance. “ROOT” represents a single derived genotype consisting of 289 SNPs with one common homozygous genotype from A. cerana (N = 7) and A. dorsata (N = 4) (32). Branches are colored to correspond to the four groups in (A). The group labeled scutellata* also includes A. m. litoria (N = 2) and A. m. capensis (N = 3). (C to I) Structure analysis (Fig. 1C, K = 4) of New World individuals. (C) Geographic distribution of Africanization in South America. Hatched region denotes hybrid zone (between latitude 31°04.993′S and 33°11.603′S). (D) Feral bees from Texas before Africanization. (E) Africanized bees (by mitochondria and morphometry) from Texas collected from 1996 to 2000. (F) Africanized bees from Arizona collected from 1996 to 2000. (G) Managed bees from Texas collected in 2005. (H and I) Time series from the Welder Wildlife Refuge (southern Texas) during Africanization (1993 to 2001). (H) Bees with African mitochondria. (I) Bees with European mitochondria.

To explore evolutionary relationships, we generated distance trees using Old World A. mellifera, and rooted these trees using an ancestral genotype derived from A. cerana and A. dorsata. The two groups that comprise European bees (M and C) were the most distantly related, despite their geographic proximity, and M was considerably more similar to A (FST = 0.242) than to either C or O (FST = 0.565 and 0.458, respectively) (Fig. 2B, table S2). Unexpectedly, the tree was rooted in Africa (Fig. 2B). Tree rooting and other relationships had strong bootstrap support and were essentially the same when we used parsimony-based methods, alternative SNP subsets, or genotypes from congeneric individuals (rather than a single derived ancestral genotype) (figs. S5 to S7). On the basis of these data, we hypothesize that A. mellifera originated in Africa and that there were at least two subsequent expansions into Eurasia—a western expansion into Europe (M) and one or more (independent) eastern expansions into Asia and Europe (O and C).

New World bees collected in Brazil (the original site of A. m. scutellata introduction in 1956) were highly Africanized, although all individuals exhibited introgression with the M group (Figs. 1C and 2C; fig. S1). Similar levels of introgression were not evident in any A. m. scutellata from Africa (Figs. 1C and 2A; fig. S1.) Two parallel transects in northern Argentina revealed transitions from predominantly African genotypes (north) to predominantly C group genotypes (south) (Fig. 2C and fig. S1B; table S2). Structure analysis revealed a consistent minority of M group alleles (Fig. 2C, black) in all individuals north and south of the hybrid zone. Retention of several M group markers has been observed at the population level (13); our results indicate that this effect occurs at the individual level. The location of the hybrid zone (Fig. 2C) appears to be stable, matching the general latitude of transition from Africanized to European bees in previous studies published in 1982 and 1991 (14, 27).

Feral bees collected in the United States before Africanization were of mixed Eurasian ancestry, dominated by C group but admixed with MandO groups (Fig. 2D and fig. S1A). After Africanization, bees defined as Africanized by morphology exhibited predominantly African ancestry admixed with M, C, and to a lesser extent, O group bees (Fig. 2, E and F; fig. S1A). During Africanization, bees from southern Texas (1993 to 2001) showed a transition from mixed European to substantially, but not exclusively, African ancestry (Fig. 2, H and I; fig. S1C). Individuals in early (1993 to 1995), middle (1996 to 1998), and late collections (1999 to 2001) with African mtDNAwere predominantly African at nuclear loci (Fig. 2H, table S2). In contrast, bees with European mtDNA became progressively differentiated from A. m. ligustica (C group) and increasingly similar to A. m. scutellata (Fig. 2I, table S2). In North America, as in South America, an unexpectedly consistent minority of M group alleles was evident in all genomes of both non-Africanized and Africanized individuals (Fig. 2, D to I, black). Recently collected bees (1999 to 2001) with African and European mtDNA were virtually indistinguishable at the SNP loci (Fig. 2, HversusI;tableS2).

To identify loci with extreme patterns of differentiation, possibly due to natural selection (28), we calculated FST among Old and New World bees at all 1136 SNP loci (Fig. 3). Pairwise differentiation between non-Africanized and Africanized bees in the New World was highly correlated with differentiation between A. m. ligustica and A. m. scutellata in the Old World (Fig. 3A) (r = 0.81), consistent with the hypothesis that Africanization involves replacement of A. m. ligustica alleles by A. m. scutellata alleles (SNPs with the highest levels of differentiation in both comparisons are indicated in Fig. 3A by filled red squares and in Table 1 and table S3 as “Set 1”). In contrast, most loci highly differentiated between A. m. scutellata and A. m. mellifera showed little differentiation between Africanized and non-Africanized bees (Fig. 3B) (r = –0.37), indicating that biased replacement of C- but not M-derived alleles [(13) and the present study] may occur throughout the genome. However, a small subset of loci that distinguish A. m. mellifera from A. m. scutellata also exhibited relatively high differentiation for Africanized versus non-Africanized New World bees (Fig. 3B, filled blue triangles). These few SNPs distinguished both C and M groups from the A group (Fig. 3C and Table 1).

Fig. 3.

FST of 1136 SNPs in introduced and native populations. (A) Differentiation in the New World (Africanized versus non-Africanized; y axis) largely corresponds to differentiation between A. m. scutellata and A. m. ligustica (x axis). Filled red squares represent the 1% of SNPs exhibiting the highest differentiation in both axes (Table 1 and table S3, “Set 1”). (B) SNPs exhibiting high differentiation between A. m. scutellata and A. m. mellifera (x axis) exhibited little differentiation in the New World (y axis). Filled blue triangles represent the 1% of SNPs exhibiting the highest differentiation in both axes (Table 1 and table S3, “Set 2”). (C) Few SNPs exhibit high differentiation in both C and M lineages relative to A, consistent with independent derivation of eastern and western European subspecies from Africa. Set 2 loci generally exhibit high differentiation in the New World (B, filled blue triangles) and in the two Old World lineages independently derived from Africa (C, open blue triangles). Several SNPs occurring in or near genes are indicated (abbreviations as in Table 1). “Africanized” refers to all individuals in Fig. 2, C (north of hybrid zone), E, and F. “Non-Africanized” refers to all individuals in Fig. 2, C (south of hybrid zone) and D.

Table 1.

Candidate loci for selection. Subset of outlier SNPs from Fig. 3A (Set 1) and Fig. 3B (Set 2) (see table S3). The frequency of predominant A. m. scutellata allele is indicated. Syn, synonymous codon change. Groups for South and North America are as indicated in Fig. 2 (arranged by increasing Africanization from left to right). UTR, untranslated region; UV, ultraviolet; GTPase, guanosine triphosphatase; PKG, cGMP-dependent protein kinase.

Old WorldSouth AmericaNorth America
SNP ligusica mellifera scutellata caucasica SouthHybrid zoneNorthEuropean1993-19951996-19981999-2001AfricanizedNear or affected geneSNP position or effectGene product
Set 1
est6550 0.00 0.95 1.00 1.00 0.06 0.83 0.99 0.21 0.38 0.64 0.88 0.91 GB11704 5′ UTR
ahb7495 0.00 0.94 0.97 0.00 0.13 0.42 0.97 0.08 0.29 0.50 0.46 0.91 GB10583 (nAchRα3) View inline Nicotinic acetylcholine receptor, α3 subunit
est8764 0.00 1.00 0.97 1.00 0.25 0.75 1.00 0.13 0.59 0.69 0.73 0.86 GB10830 S→A
est9209View inline 0.03 0.98 1.00 0.36 0.19 0.75 0.97 0.21 0.54 0.62 0.79 0.91 GB10514(Tubα1) Syn Tubulin, α1 chain
est9211View inline 0.03 0.98 0.97 0.32 0.19 0.75 0.97 0.21 0.55 0.67 0.79 0.95 GB10514 (Tubα1) Syn Tubulin, α1 chain
Set 2
ahb12140 0.00 0.00 1.00 0.00 0.00 0.33 0.76 0.04 0.20 0.38 0.44 0.73 GB18171 (UVop) Syn UV-sensitive opsin
ahb11258 0.00 0.00 0.92 0.00 0.00 0.25 0.74 0.00 0.21 0.40 0.65 0.86 GB15150 Intron GTPase activator
ahb11774 0.08 0.03 0.94 0.00 0.00 0.50 0.76 0.04 0.23 0.55 0.54 0.50 GB15050 ∼1500 bp 5′
est424 0.25 0.20 0.97 0.00 0.19 0.92 0.92 0.17 0.41 0.45 0.71 0.86 GB19539 Syn
est10185 0.00 0.00 0.91 1.00 0.06 0.50 0.81 0.13 0.32 0.57 0.58 0.45 GB18394 (for) 3′ UTR foraging; PKG
ahb9731 0.00 0.05 0.78 0.00 0.00 0.17 0.69 0.00 0.18 0.57 0.50 0.77 GB15214 Syn
  • View inline* These two SNPs are 939 base pairs (bp) apart.

  • View inline Approximately 3000 bp 5′ of start codon for gene prediction GB10583; this gene has an alternate predicted structure (S.C_Group5.24000030B) (3) that places the SNP 30 bp upstream of the start codon.

  • Of 19 SNPs identified as potential sites for selection (Fig. 3), 11 occurred in or near genes (Table 1 and table S3). Six occurred in coding regions, but only one was predicted to cause a nonsynonymous amino acid change (in a gene of unknown function; Table 1). Although some of these SNPs may alter gene function, it is also possible that some or all have hitchhiked with linked polymorphisms that are under selection.

    These data provide a global view of the population genetic structure of A. mellifera. Future studies can capitalize on the large number of SNPs presented here to further investigate population genetic patterns within and among populations and genes under selection.

    Supporting Online Material

    Materials and Methods

    Tables S1 to S3

    Figs. S1 to S8


    References and Notes

    View Abstract

    Stay Connected to Science

    Navigate This Article