The expansion of Africanized honeybees from South America to the southwestern United States in <50 years is considered one of the most spectacular biological invasions yet documented. In the American tropics, it has been shown that during their expansion Africanized honeybees have low levels of introgressed alleles from resident European populations. In the United States, it has been speculated, but not shown, that Africanized honeybees would hybridize extensively with European honeybees. Here we report a continuous 11-year study investigating temporal changes in the genetic structure of a feral population from the southern United States undergoing Africanization. Our microsatellite data showed that (1) the process of Africanization involved both maternal and paternal bidirectional gene flow between European and Africanized honeybees and (2) the panmitic European population was replaced by panmitic mixtures of A. m. scutellata and European genes within 5 years after Africanization. The post-Africanization gene pool (1998–2001) was composed of a diverse array of recombinant classes with a substantial European genetic contribution (mean 25–37%). Therefore, the resulting feral honeybee population of south Texas was best viewed as a hybrid swarm.
THE evolutionary significance of natural hybridization has been debated for decades (Mayr 1942; Anderson 1949; Harrison 1993; Arnold 1997). At one extreme, it has been argued that natural hybridization is an evolutionary dead end due to formation of inviable and/or infertile hybrids (Mayr 1942; Barton and Hewitt 1985, 1989). At the other extreme, it has been suggested that natural hybridization may lead to new evolutionary lineages due to formation of relatively fit hybrids that expand into novel habitats (Anderson 1948; Arnold 1997; Ellstrand and Schierenbeck 2000; Bleeker 2003). A third potential evolutionary outcome is expansion of an intermixed form within the resident progenitor's habitat, in which case the degree of mixing between hybridizing forms may range from formation of a hybrid swarm to genetic assimilation of one form by the other (Childs et al. 1996; Rhymer and Simberloff 1996; Perry et al. 2001). Natural hybridization and introgression have been reported in a growing number of biological invasions (Rhymer and Simberloff 1996). Such invasive events are of great interest to the discipline of evolutionary biology because they provide unique opportunities to study evolutionary processes at initial stages of secondary contact of divergent genomes.
Our study of hybridization, which deals with one of the most spectacular biological invasions yet documented, that of the Africanized honeybees (Apis mellifera L.), builds on the previous studies of Africanization. Africanized honeybees derive from a founder population of the tropical-evolved African subspecies Apis mellifera scutellata brought from South Africa into Brazil in 1956 to interbreed with previously imported temperate-evolved European honeybees. The purpose of the introduction was to create a honeybee better adapted to tropical conditions, because the European honeybee colonies were poor honey producers and would survive only when intensively managed (Nogueira-Neto 1964; Kerr 1967). Soon after the introduction, queens of A. m. scutellata were accidentally released into the natural environment (Spivak et al. 1991) and their descendants have since expanded throughout South and Central America and have established large feral populations where European colonies could not thrive (Taylor 1988; Roubik and Boreham 1990; Winston 1992). The leading edge of the expanding front reached the southern United States in 1990, only 33 years after initial release (Sugden and Williams 1990). Prior to arrival of Africanized honeybees, the United States sustained large feral and managed populations of honeybees predominantly derived from eastern (A. m. ligustica, A. m. carnica, and A. m. caucasia) and western (A. m. mellifera) European subspecies (Rubink et al. 1990; Schiff et al. 1994; Schiff and Sheppard 1995, 1996; Loper et al. 1999). With the arrival and subsequent expansion of Africanized honeybees into Texas, New Mexico, Arizona, Nevada, and California, the European honeybee gene pool in those areas is being converted to one with a substantial African component (Rubink et al. 1996; Loper et al. 1999; Pinto et al. 2004).
The underlying mechanism for dispersal and the genetic composition of Africanized honeybee populations has been intensively studied in the American tropics (Del Lama et al. 1988; Hall and Muralidharan 1989; Lobo et al. 1989; Smith et al. 1989; Hall 1990; Rinderer et al. 1991; Lobo and Krieger 1992; Quezada-Euán 2000; Clarke et al. 2002) and to a lesser degree in more temperate regions of South (Sheppard et al. 1991a; Diniz et al. 2003; Quezada-Euán et al. 2003) and North America (Loper et al. 1999; Pinto et al. 2004). Studies of neotropical feral populations suggest that Africanized honeybees expanded by maternal migration (Taylor 1988; Hall and Muralidharan 1989; Smith et al. 1989; Hall and Smith 1991) and have introgressed unequal nuclear and mitochondrial proportions of European genes (Hall 1990; Sheppard et al. 1991b; Hall and McMichael 2001; Clarke et al. 2002). Surveys of allozymes (Del Lama et al. 1988; Lobo et al. 1989; Sheppard et al. 1991b; Diniz et al. 2003), restriction fragment length polymorphisms (Hall 1990; McMichael and Hall 1996; Hall and McMichael 2001), randomly amplified polymorphic DNA (Suazo et al. 1998), and microsatellites (Clarke et al. 2002) showed low-to-intermediate levels of European alleles in the Africanized gene pool, whereas surveys of mitochondrial DNA (mtDNA) showed fixed-to-low levels of African-type mitochondria (Hall and Muralidharan 1989; Smith et al. 1989; Hall and Smith 1991; Sheppard et al. 1991b; Clarke et al. 2001; Diniz et al. 2003). The varying proportions of European markers found in different studies seem to be dependent on (1) the size of the preexisting European honeybee population; (2) the time since secondary contact; and possibly (3) the nature of the population, i.e., managed vs. feral (Hall and Muralidharan 1989; Smith et al. 1989; Rinderer et al. 1991; Clarke et al. 2002). The overall pattern of Africanization that emerged from these studies is that Africanized honeybees may hybridize extensively when (managed) European honeybees are common. However, over time European alleles largely disappear mainly because of continued migration of African genotypes into the area and possibly because hybrid workers exhibit reduced fitness (Harrison and Hall 1993; Schneider and Hall 1997; Schneider et al. 2003).
As the Africanized honeybees are expanding in southwestern United States, a hybrid zone, which potentially mirrors a hybrid zone of South America, is forming. Studies of the stable hybrid zone in South America reported increasing and persisting proportions of European genes in north-south clines. In the subtropical-temperate transitional areas, a diverse array of associations of European and African markers was found, suggesting that the breakdown of the African genome was occurring (Lobo et al. 1989; Sheppard et al. 1991a; Diniz et al. 2003). While the hybrid zone in North America has not been studied, it is possible that its dynamics have been changed by the parasitic mite Varroa destructor. Since introduction in 1987 (De Guzman et al. 1997), Varroa has caused severe losses in managed and feral European populations (Kraus and Page 1995; Loper et al. 1999; Pinto et al. 2004). Further, several studies have shown that European honeybees are more susceptible to Varroa than Africanized honeybees (Moretto et al. 1991, 1993; Message and Gonçalves 1995; Guzman-Novoa et al. 1996). Therefore, the combined effect of a depauperate population and differential susceptibility to Varroa could have changed the rate and extent of hybridization between European and Africanized honeybees.
We completed an 11-year study of mtDNA of a feral population of European honeybees from the southern United States (Texas) undergoing successive invasions of both the Africanized honeybee and Varroa. This population displayed a rapid decline in the number of colonies of honeybees possessing non-A. m. scutellata mitotypes, but not in colonies of honeybees with the A. m. scutellata mitotype (Pinto et al. 2004). In spite of the initial trend, the authors found a substantial contribution of non-A. m. scutellata mitotypes several years after both invasions, suggesting a persistence of European genes as predicted in a more temperate setting. While the finding suggested that the feral population lies within a hybrid zone, it did not rule out complete genetic assimilation of the European genome by Africanized honeybees. In the present study, nuclear DNA (microsatellites) markers were screened from the 11-year honeybee collection and used concurrently with the mtDNA to better understand the underlying mechanisms of Africanization and the genetic composition of the feral population over time. This study reports, on a fine temporal scale, the pattern of hybridization and introgression involved in this event by examining estimates of single-locus disequilibria, linkage disequilibria, cytonuclear associations, and admixture proportions. This is the first comprehensive study concurrently using nuclear and mtDNA markers of a feral population in a more temperate setting encompassing the whole invasion process.
MATERIALS AND METHODS
A total of 335 honeybee workers was sampled from feral colonies between 1991 and 2001 at the Welder Wildlife Refuge (WWR), which is located at 28° N latitude, San Patricio County, Texas. Each honeybee worker represented a different and a new colony; i.e., no colonies were resampled. The 11-year collection comprised honeybees from colonies obtained from tree cavities (244 samples) and swarm traps (91 samples) (see Pinto et al. 2004 for details on sampling procedure).
Honeybee workers of two additional populations (one European derived and one A. m. scutellata derived) were sampled for use as reference populations. The European-derived population (referred to as “Eur.” throughout this article) was sampled from 50 feral colonies caught in swarm traps between 1988 and early 1990 (prior to Africanization). The swarm traps were deployed across an east-west 120-km-long transect, established in the lower Rio Grande Valley (between the cities of Brownsville and La Joya) of South Texas (Rubink et al. 1996). The A. m. scutellata-derived population from Brazil (referred to as “Braz.” throughout this article) was sampled in 2002 from 43 colonies maintained in the apiary of the University of São Paulo (Ribeirão Preto, São Paulo, Brazil). The sampled colonies, for the most part, originated as feral swarms collected in the Luiz Antonio area (∼100 km north of the initial release site of A. m. scutellata in Rio Claro). The swarms were maintained without subsequent queen management practices and therefore were considered representative of the local feral population (D. De Jong, personal communication). We assumed that the sample from South Texas was representative of the pre-Africanized WWR population. The sample from Brazil was our best representative of the A. m. scutellata-derived founder population and was regarded as the source of the migrants. The honeybees were either frozen on dry ice (samples from South Texas and WWR) or stored in 95% ethanol (samples from Brazil) and transported back to the laboratory where they were stored at −80° and 4°, respectively, until DNA extraction.
Total DNA was extracted from the honeybee mesosoma (thorax plus the first fused abdominal segment) using a QIAamp DNA mini kit (QIAGEN, Valencia, CA) and stored at −20° until analyzed. A total of 428 honeybee samples were scored at 12 microsatellite loci (A14, A7, A88, A107, A113, A35, A28, A79, A43, A8, IM, and ED1; see Estoup et al. 1994, 1995; Oldroyd et al. 1995; and Rowe et al. 1997 for primer sequences). Single polymerase chain reaction (PCR) amplifications were performed in 8 μl total volume containing 0.5× Taq DNA polymerase buffer, 1.5 mm MgCl2, 0.2 mm of each dNTPs, 2 pm of each primer, 1 μl of template DNA, and 0.4 units of Taq DNA polymerase (Promega, Madison, WI) for all primer pairs. The forward primer for each marker was labeled with one of three fluorescent dyes (HEX, 6-FAM, or NED; Applied Biosystems, Foster City, CA). The PCR conditions consisted of an initial 3-min denaturation step at 94° followed by 30 cycles of 94° for 15 sec, annealing temperature at 55° (except for A79 at 60°) for 15 sec, and 72° for 5 sec, plus a final extension at 72° for 10 min. The PCR products were loaded on Long Ranger polyacrylamide gels (BioWhittaker Molecular Applications, Rockland, ME) and run on an ABI Prism 377 automated DNA sequencer using HD Rox 400 (Applied Biosystems) as the internal size standard. The software Genotyper 2.5 (Applied Biosystems) was used for allele identification and comparison.
Genetic diversity was assessed using unbiased estimates of gene diversity (Nei 1987) and allelic richness (Petit et al. 1998) for each microsatellite locus in each sample. Allele frequencies, number of alleles (A), proportion of heterozygotes (Hp), unbiased gene diversity (Hd), and allelic richness (Rs, a measure of the number of alleles independent of sample size) per locus and population sample were computed using the FSTAT version 2.9.3. package (Goudet 2001). Differences in average unbiased gene diversity and allelic richness between pairs of samples were assessed by Wilcoxon's signed rank test (Snedecor and Cochran 1978) using SPSS version 11.0 (Norušis 1993).
Departure from Hardy-Weinberg equilibrium (HWE) was tested for each locus and population sample using the Hardy-Weinberg exact test. A global probability value over all loci was obtained for each population sample using the procedure of Raymond and Rousset (1995a). Heterozygote excess and deficiency were tested at each locus and across loci for each sample using the score test (U-test) according to Rousset and Raymond (1995). Fisher's exact tests were performed to assess genotypic linkage disequilibrium for all locus pairs within each population sample (Raymond and Rousset 1995a). Population sample pairwise comparisons were performed to test for genic (Raymond and Rousset 1995b) and genotypic (Goudet et al. 1996) differentiation across loci. These analyses were performed by the Markov chain method (Guo and Thompson 1992) with 10,000 dememorization steps, 1000 batches, and 10,000 iterations per batch using GENEPOP version 3.4 (Raymond and Rousset 1995a).
Reference and WWR population sample pairwise values for FST were computed according to Weir and Cockerham (1984) using GENEPOP and then departures from zero were tested by 1000 permutations using GENETIX version 4.04 (Belkhir et al. 2002). Distance measures, termed proportion of shared alleles (DPS, Bowcock et al. 1994), between pairs of populations were generated using the program MSAT. These distances were used to build a neighbor-joining tree (Saitou and Nei 1987) using PAUP version 4.0b 10 (Swofford 2002). Where applicable throughout the analysis, statistical significance levels were adjusted for multiple comparisons using the sequential Bonferroni procedure to correct for type I error (Rice 1989).
The European and Brazilian reference population samples were used as parentals to compute admixture proportions in the WWR population using the estimators mR and mY. The least-squares estimator mR is based on a comparison of allele frequencies between parental and admixed populations and assumes identity by descent (Roberts and Hiorns 1965). The estimator mY is based on a coalescent approach that takes into account molecular information as well as allele frequencies. This estimator, which assumes the single stepwise mutation model, was computed using the squared difference in allele size as molecular distance. Both estimators were computed using the program ADMIX version 1.0. The bootstrap average and bootstrap standard deviation of the admixture coefficients were estimated over 1000 replications (Bertorelle and Excoffier 1998). Differences in mR and mY over time were assessed by a paired t-test using SPSS.
The European and Brazilian reference samples were also used as parentals to simulate 100 sets of F1, European backcross, and Brazilian backcross population samples of 50 individuals each. The multilocus genotype of each simulated individual was generated by randomly taking one allele per locus from each of the parental populations, according to their frequencies, to form the respective “gametes” using a program written by J. B. Patton (personal communication). The simulated F1 and reciprocal backcrosses were used as additional reference samples to build the above-mentioned neighbor-joining tree and to classify individual honeybees from WWR over time using GENECLASS (Cornuet et al. 1999), a software package that assigns individuals of unknown origin to a set of reference samples using either likelihood (frequency or Bayesian) or distance-based methods. To compare the performance of frequency, Bayesian, Cavalli-Sforza and Edwards (1967) chord distance, and (δμ)2 Goldstein et al. (1995) distance methods, multilocus genotypes of 1000 European, 1000 Brazilian, 1000 F1, 1000 European backcross, and 1000 Brazilian backcross individuals were randomly generated by crossing European × European, Brazilian × Brazilian, European × Brazilian, F1 × European, and F1 × Brazilian, respectively. Using one random set (among 100) of reference samples, each simulated individual was assigned and the percentage of misclassification was calculated for each method. The most accurate method was examined over all 100 sets of simulated reference samples and then used to classify individual honeybees from the WWR. Assignment scores for each WWR honeybee were computed 100 times, using the 100 sets of simulated reference samples. Classification of individual honeybees was based on the predominant assignment given for the 100 replications.
Allele frequencies, number of alleles, allelic richness, proportion of heterozygotes, and unbiased gene diversity are given for each locus and population sample in supplementary data at http://www.genetics.org/supplemental/. All loci showed high variability with the total number of alleles ranging from 11 (A8) to 38 (A7), the mean number of alleles ranging from 5.8 (WWR 1991) to 15.1 (WWR 1995), the mean allelic richness ranging from 5.833 (WWR 1991) to 9.168 (WWR 2000), the mean proportion of heterozygotes ranging from 0.700 (WWR 1994) to 0.876 (WWR 1999), and finally the mean gene diversity ranging from 0.681 (WWR 1991) to 0.874 (Brazilian) across all population samples.
Genetic diversity was significantly higher in the Brazilian reference sample (Rs = 8.890; Hd = 0.874) than in the European (Rs = 5.934; Hd = 0.693) reference sample (P = 0.002 for both allelic richness and gene diversity). Pairwise comparisons between the reference and WWR samples showed, first, that mean allelic richness and mean gene diversity were significantly lower in the European (Rs = 5.934; Hd = 0.693) than in the WWR after 1994 (Rs = 6.653–9.168) and 1995 (Hd = 0.800–0.874), respectively (P ≤ 0.005); and, second, that mean allelic richness and mean gene diversity were significantly higher in the Brazilian (Rs = 8.890; Hd = 0.874) than in the WWR (P ≤ 0.006) prior to 1994 (Rs = 5.833–6.653) and 1995 (Hd = 0.681–0.800), respectively. No significant differences in genetic diversity were found for the remaining reference-WWR sample pairwise comparisons (P ≥ 0.060). During the 11-year study, the WWR population experienced a striking increase in gene diversity and allelic richness, with the greatest increment observed between 1994 and 1996 (Figure 1).
Hardy-Weinberg equilibrium and genotypic linkage equilibrium:
FIS estimates and respective significance, P-values for HWE, P-values for heterozygote deficiency, and P-values for heterozygote excess and their standard errors per locus are shown in supplementary data at http://www.genetics.org/supplemental/ and per population sample in Table 1. In the reference European sample, no significant departures from HWE were detected for any locus and across loci (0.085 ≤ P ≤ 0.963). Individual and multilocus U-tests for heterozygote deficiency and heterozygote excess yielded a single significant value at A107 (P = 0.036 for heterozygote excess), which became nonsignificant following sequential Bonferroni correction. In contrast, the Brazilian reference sample showed a deviation from HWE (P = 0.013), caused by a deficiency of heterozygotes (P = 0.001) across all loci. Single-locus P-values obtained with Fisher's exact test, U-tests, and/or FIS revealed a deficiency of heterozygotes at A88, A14, A35, and A79 and an excess of heterozygotes at ED1 in this reference population sample (supplementary data at http://www.genetics.org/supplemental/). However, when applying sequential Bonferroni correction, only A35 and A79 exhibited a significant deficiency of heterozygotes (P = 0.000). Fisher's exact tests for genotypic linkage disequilibrium between loci pairs produced eight (four for the European and four for the Brazilian reference sample) significant P-values from 132 comparisons. None of the loci pairs were shared by the two populations and none of the P-values was significant following sequential Bonferroni's correction.
In the WWR population, significant departures from HWE across all loci were detected in 1995, 1996, and 2000 (0.002 ≤ P ≤ 0.005; Table 1). These deviations were caused by a deficiency of heterozygotes (0.0001 ≤ P ≤ 0.007). In addition to those years, multilocus U-tests revealed a deficiency of heterozygotes in 1994, 1997, and 2001 (0.000 ≤ P ≤ 0.041) and an excess of heterozygotes in 1991 (P = 0.029). However, when applying the sequential Bonferroni procedure, neither 1991 nor 2001 WWR population samples deviated from HWE. At the single-locus level, departures from HWE occurred at all loci except ED1 for at least one of the sampling years, as indicated by Fisher's exact test, U-tests, and/or FIS P-values (supplementarydataathttp://www.genetics.org/supplemental). The greatest number of significant P-values was observed in 1995 (8 values of P < 0.05 given by Fisher's exact test, U-tests, and FIS) and 1996 (8 values of P < 0.05 given by Fisher's exact test, U-tests, and FIS). The majority of the 16 P-values were marginally significant and subsequently became nonsignificant when corrected for type I error. Across all years, locus A79 exhibited the most consistent deficiency of heterozygotes, which was detected by Fisher's exact test, U-test, and/or FIS in 1995, 1996, 1997, 1999, 2000, and 2001. The hypothesis of null alleles at this locus was dismissed because no heterozygote deficiency was detected in the reference European and early WWR population samples. However, it was detected in the reference Brazilian. Because deficiency of heterozygotes seems to be correlated with Africanization, this locus deserves further investigation. Fisher's exact tests for linkage disequilibrium yielded 46 significant values from 704 comparisons (35 significant values are expected at the 5% level). The most significant values were in the 1995 (22) and 1994 (9) population samples (Figure 2). When applying sequential Bonferroni correction, only 5 (1 in 1994 and 4 in 1995) of the 46 above-mentioned P-values were significant. Two pairs of loci, A28/ED1 (mapping distance 23 cM) and A7/A14 (mapping distance 22 cM), are genetically linked, while the eight remaining loci are unlinked (M. Solignac, personal communication). However, all 12 loci behaved similarly and no single-locus pair produced consistently significant values for all years after population mixing. Across the 11-year study, the mean P-value was lowest in 1994 and 1995. P-values obtained from 1996 onward suggested a rapid dissipation of linkage disequilibrium (Figure 2).
The multilocus distance measure (DPS), multilocus FST, and tests for homogeneity of genic and genotypic distributions (genic and genotypic differentiation) showed a high level of differentiation between Brazilian and European reference samples (P < 0.000 for genic and genotypic differentiation and FST; Tables 2 and 3). Pairwise comparisons between reference European and WWR samples revealed an increase in genic and genotypic differentiation over time, with FST values becoming significantly different from zero since 1995. In contrast, pairwise comparisons between reference Brazilian and WWR population samples showed a decrease in differentiation over time. In spite of the trend, all pairwise comparisons yielded significant values, suggesting that the WWR population was genetically distinct from the Brazilian. Values of DPS, FST, and genic and genotypic differentiation indicated that the WWR population samples could be divided into pre- and post-1996 populations. Between 1991 and 1994 pairwise comparisons of genic and genotypic differentiation were nonsignificant, and pairwise multilocus FST values were not significantly different from zero. The same pattern was shown by pairwise comparisons performed between 1997 and 2001 samples. In contrast, pre-1996 vs. post-1996 sample comparisons showed a high level of differentiation, suggesting that a profound genetic change occurred in a short time frame. A neighbor-joining tree confirmed this pattern (Figure 3). Until 1996, the WWR samples showed closer genetic affinities with the reference European sample. After 1996, the WWR population samples were clustered closer to the reference Brazilian population.
Both mR and mY estimators revealed a striking change in the admixture proportions of the WWR population over time (Figure 4). For all years, mY estimated a greater contribution of the parental Brazilian to the WWR population than did mR (P = 0.000). The proportion of A. m. scutellata introgressed alleles in the WWR population increased from 2 to 11% in 1991 to 63–68% in 2001 with mR and mY estimators, respectively. The steepest increase in the proportion of introgressed A. m. scutellata alleles occurred between 1994 (mR = 9.6%; mY = 12.1%) and 1997 (mR = 52.5%; mY = 63.4%). From 1998 on, the nuclear proportion of A. m. scutellata genes in the WWR population appeared to have stabilized at 62.8% ± 0.042 (SD) for mR and 74.7 ± 0.043 (SD) for mY. Interestingly, as shown in Figure 4, the temporal pattern of introgression of nuclear genes was paralleled by A. m. scutellata mitotype (Pinto et al. 2004).
Associations of nuclear and mtDNA markers:
Among the assignment methods that were tested in GENECLASS, Bayesian showed the best performance with the highest rate (68.3%) of individuals correctly assigned (Table 4). Cornuet et al. (1999) also found this method to be the most accurate even when the assumptions of HWE and genotypic linkage equilibrium were violated. The percentage of misclassification given by the Bayesian method (and the other methods tested; data not shown) varied across types of simulated individuals. More than 96% of the simulated European and Brazilian individuals were correctly assigned whereas <50% of backcross individuals were correctly assigned (Figure 5). The average assignment scores (values produced by the Bayesian methodology of Rannala and Mountain 1997 as modified in GENECLASS), which were obtained for simulated individuals, are shown in Figure 6. The highest assignment scores were produced by reciprocal backcross (19.04 ± 2.77 for European backcross; 22.15 ± 2.39 for Brazilian backcross) and F1 (21.11 ± 2.45) individuals whereas the lowest were produced by European (13.64 ± 1.75) and Brazilian (20.03 ± 1.49) individuals. As crosses became more complex, the overlap among scores increased and the power of the assignment test decreased (Figure 6). The average assignment score obtained for the WWR population increased steeply between 1994 and 1996 (Figure 7), which suggests an increase in hybrid complexity of multilocus genotypes. After 1997, the average assignment score stabilized at high values. Comparisons of ranges of the assignment scores for simulated (Figure 6) and WWR individuals (Figure 7) demonstrate that crosses were characteristically more complex than simple backcrosses.
Figure 8 shows the temporal frequency of non-A. m. scutellata (eastern European, western European, and A. m. lamarckii; Pinto et al. 2004) and A. m. scutellata mitotypes and corresponding nuclear assignment classes (European, Brazilian, F1, European backcross, or Brazilian backcross). The proportion of individuals of eastern European, western European, A. m. lamarckii, and A. m. scutellata mitotypes that were genotypically closer to European, Brazilian, F1, European backcross, or Brazilian backcross reference population changed radically in WWR over time. Prior to 1994 most individuals exhibited a non-A. m. scutellata mitotype and were classified as either European or (less frequently) European backcrosses. Between 1994 and 1996, individuals of non-A. m. scutellata maternal descent were genotypically more homogeneous than those with A. m. scutellata mitotype and an association between nuclear and mtDNA genes of non-A. m. scutellata origin was apparent (Figure 8). Significant heterogeneity in genic and genotypic distributions between both groups of mitotypes was found in that 3-year period (0.000 ≤ P ≤ 0.005; the three non-A. m. scutellata mitotypes were pooled to increase sample size for added statistical power). However, a breakdown of non-A. m. scutellata cytonuclear associations appeared to have occurred as no significant differences in genic and genotypic distributions were found from 1997 onward (0.117 ≤ P ≤ 0.932; the three non-A. m. scutellata mitotypes were pooled). After 1997, most individuals exhibited a mixed ancestry with no observed association between mitotype and nuclear composition.
While Figure 8 allows a detailed look at mitotype by nuclear genome interactions, it does not allow a concise summary statistic of those interactions. Therefore, we generated Figure 9 using the multilocus FST values of the genomic composition of the WWR honeybees characterized by possession of either the A. m. scutellata or the non-A. m. scutellata mitotype by year as compared to the two reference populations. The nuclear genomes of individuals possessing either the A. m. scutellata or the non-A. m. scutellata mitotype began to merge at least as early as 1995. By 1997 the nuclear genomic composition associated with the respective mitotypes had converged to a point where the multilocus test for both genic and genotypic distributions was nonsignificant as previously discussed. Multilocus FST values for comparisons of the years from 1998 to 2001 virtually overlay one another, indicating the nonexistence of cytonuclear associations (Figure 9).
Mechanism of Africanization:
mtDNA data revealed the year 1993 as the time of Africanization onset in the WWR population (Pinto et al. 2004), whereas estimates of nuclear racial admixture indicated a contribution of A. m. scutellata alleles beginning in 1991 (Figure 4). This result suggests that Africanization could have started prior to 1993 through matings of migrant drones with resident European queens. Alternatively, the early evidence of admixture could be an artifact of the analysis because for 1991 the error bars for the estimated admixture proportion, as calculated by both mR and mY, included zero. Hunter et al. (1993) documented the presence of morphologically defined Africanized honeybees in a county adjacent to San Patricio county, where WWR is located, in 1991. Therefore, it is possible that Africanized drones may have hybridized with European queens, which established colonies at WWR as early as 1991. In 1992 the error bars around the mR estimator do not overlap with zero, whereas those of the mY estimator do. While the results are not entirely statistically significant, the observation of 14 swarms of Africanized honeybees in San Patricio county in 1992 (Hunter et al. 1993) supports the contention that paternal gene flow did occur and was the likely cause for the African component of the admixture proportions. Additional support is given by Figure 8, which shows an excess of individuals assigned as European backcrosses for both 1991 and 1992 relative to the simulations to estimate the frequency of misassignment illustrated in Figure 5 for the European reference population. To conclude, the body of evidence suggests that Africanization in WWR started by matings of Africanized drones rather than by immigration of Africanized queens.
In the neotropics, genetic studies show that in feral populations maternal gene flow was the primary driving force of Africanization (Taylor 1988; Hall and Muralidharan 1989; Smith et al. 1989; Hall and Smith 1991), whereas in managed populations the primary force was paternal gene flow (Rinderer et al. 1991; Clarke et al. 2002). In contrast with the neotropics, where a virtually “pure African” expanding front was observed, migrants arriving at WWR were of mixed ancestry (Figure 8). This had the effect of diluting the input of the African nuclear loci relative to the African mtDNA, which could account for maintenance of Hardy-Weinberg equilibrium prior to 1994. Regardless of the mechanism of initial Africanization, the most dramatic genetic change of the WWR population occurred after the arrival of the Africanized queen in 1993.
Turnover in the population:
A 4-year period of Hardy-Weinberg disequilibrium due to deficiency of heterozygotes was observed in WWR between 1994 and 1997 (Table 1). Hardy-Weinberg disequilibrium was anticipated upon arrival of migrants as mixing of distinct gene pools prior to reproduction would cause a Wahlund effect. While the Wahlund effect is expected to disappear after one generation of random mating, the persistence of heterozygote deficiency for 4 years could have been due to: (1) positive assortative mating, (2) continuing influx of genetically distinct migrants, and/or (3) selection. In early reports of Africanization, positive assortative mating, due to asynchrony in mating flight times, was suggested as an explanation for apparent partial reproductive isolation between European and Africanized honeybees (reviewed in Michener 1975). This hypothesis has long been disregarded because findings of associations between European and African mitotypes with a range of morphological and allozymic phenotypes demonstrated interbreeding between the honeybee types (Rinderer et al. 1991; Sheppard et al. 1991a). The patterns depicted in Figures 8 and 9 clearly show bidirectional gene flow in the WWR population. As a result of gene exchange, associations between non-A. m. scutellata mitotypes and European nuclear classification gradually decayed and by 1997 there were no nuclear multilocus genotypic and genic differences between non-A. m. scutellata and A. m. scutellata mitotypes (P ≥ 0.117).
The A. m. scutellata mitotype increased in frequency between 1993 and 1995 from 4.5 to 28.6% (Figure 4). This could have been due to in situ reproduction of recently arrived migrants and immigration of additional Africanized swarms that further contributed to heterozygote deficiency. By 1997 the genetic composition of the population had become predominantly A. m. scutellata. The turnover resulted primarily from a dramatic reduction in the number of colonies of non-A. m. scutellata matrilines that was not observed for colonies of A. m. scutellata matrilines (Pinto et al. 2004). The authors attributed this differential loss to the parasitic mite V. destructor, which was first detected in WWR in early 1995. The observed mortality was consistent with the observation that European honeybees are more susceptible to Varroa than Africanized honeybees (Moretto et al. 1991, 1993; Message and Gonçalves 1995; Guzman-Novoa et al. 1996). In a crossbreeding experiment between Africanized and European honeybees Guzman-Novoa et al. (1996) found that (1) susceptibility to becoming infested by Varroa was least likely in the brood of Africanized honeybees, followed by European, F1 of European mother, and F1 of Africanized mother, respectively, and (2) adult European honeybees were more likely to become infested with Varroa than were adult Africanized and hybrid honeybees. In conclusion, because of the differential susceptibility of honeybees to Varroa and because the WWR population was hit at an early stage of Africanization, when European and likely F1 individuals were abundant, we suggest that differential selection pressures, possibly combined with the ongoing influx of migrants and some assortative mating, contributed to the observed disequilibrium throughout 1996 and 1997.
The temporal pattern of linkage disequilibrium experienced by the WWR population, as given by the average P-values of Fisher's exact test, is shown in Figure 2. Contrary to expectations, linkage disequilibrium decayed more rapidly than heterozygote deficiency. This circumstance is even more striking, given that two pairs of loci are genetically linked (A28/ED1 and A7/A14; M. Solignac, personal communication). Typically, more generations are required to reach linkage equilibrium than to reach Hardy-Weinberg equilibrium. In addition, both equilibria are approached at a slower pace when mating is assortative (Hartl and Clark 1997).
It is possible that the rapid dissipation of linkage disequilibrium (LD) results from Fisher's exact test for LD being less powerful than the U-test for heterozygote deficiency in detecting small deviations from expectations. Reduced statistical power for the Fisher's exact tests is not unexpected, given that European and A. m. scutellata subspecies share relatively common alleles at several loci (supplementary data at http://www.genetics.org/supplemental/; Estoup et al. 1995; Clarke et al. 2002) and given the moderate sample size (especially from 1996 onward), together with high levels of polymorphism. Nonetheless, given that Hardy-Weinberg proportions were reached in 1998 and linkage disequilibrium decayed even more quickly (Figure 2), heterozygotes (intermating between migrants and resident honeybees) must have been formed in significant numbers.
Mechanism of population turnover:
In the American tropics, feral colonies of Africanized honeybees that exhibited little or no introgression of European genes are reported to have built up to high densities within 2 to 3 years after initial colonization (Taylor 1985). It has been suggested that rapid expansion of a nearly “pure” neotropical African population has resulted from (1) high fitness of Africanized honeybees exhibiting virtually no introgression of European alleles, (2) reduced competition from a small and poorly adapted European population, and (3) cytonuclear incompatibility in hybrids of European mothers (Taylor 1985, 1988; Rinderer et al. 1991; Sheppard et al. 1991a; Harrison and Hall 1993; Schneider et al. 2003).
The WWR, which lies near the subtropical-temperate interface, supported a high density of European colonies at the arrival of the first wave of Africanized swarms (11.7 colonies/km2 in 1995; Baum 2003). Therefore, our prediction was that rate and extent of Africanization in the WWR population would be slower and lower than in neotropical populations. We expected (1) increased fitness of European honeybees under these more temperate climatical and ecological conditions, (2) higher frequency of matings with European honeybees due to the larger resident population size, and (3) strong competition of feral European colonies with the Africanized colonies for the density-limited food and nest resources. Five years after the initial observation of honeybees with African mitotypes at WWR, “pure” European honeybees were replaced by Africanized honeybees exhibiting varying degrees of introgressive hybridization (Figure 8). This level of hybridization approached the nuclear composition of the honeybees of Yucatan 13 years after the onset of Africanization (Rinderer et al. 1991; Quezada-Euán et al. 1996; Clarke et al. 2002). Complete and rapid replacement of resident European honeybees by Africanized honeybees in WWR would be expected if the latter exhibited higher fitness than the former. Whether European honeybees were equally or better fit than Africanized honeybees in the ecological and climatical conditions of southern Texas may never be known. Regardless of the relative fitnesses in a Varroa-free environment, Africanized and hybrid honeybees were likely more successful in the WWR environment that contained Varroa (Guzman-Novoa et al. 1996). This appears the most logical explanation for the unexpectedly steep increase in A. m. scutellata genetic proportions observed between 1996 and 1999 (Figure 4). Freed of competition for limited food and nest resources, expansion of the Africanized honeybees was greatly facilitated. During this period the available cavities were rapidly filled, approaching carrying capacity by 1999 (Baum 2003). Therefore, we suggest that the nearly coincidental arrival of Varroa with expansion of the Africanized front hastened the demise of “pure” European honeybees and thereby had a major role in restructuring the post-Africanization WWR population. Had the resident European population not collapsed, the rate of Africanization would have been slower. At the same time, hybridization may have prevented the loss of European alleles if F1 and later-generation hybrids were more tolerant of Varroa than “pure” European honeybees.
In contrast with neotropical feral populations, the nuclear and mitochondrial makeup of the WWR population changed symmetrically over time (Figure 4). The result does not support the hypothesis of cytonuclear incompatibilities in hybrids of European mothers, as was suggested to explain the paucity of European mitochondria in neotropical populations (Harrison and Hall 1993; Schneider et al. 2003). However, under the temperate conditions the selective pressure generating cytonuclear incompatibility may not be operative. Definitive conclusions regarding cytonuclear asymmetry in WWR cannot be drawn, given the very recent secondary contact. The long-term evolutionary outcome of nuclear and mitochondrial introgression is uncertain because European alleles may be disfavored in an A. m. scutellata genetic background. If negative selection is occurring, then European alleles may gradually be eliminated from the WWR population. Alternatively, if selection favors individuals of mixed ancestry, then the European contribution will persist.
We suggest that the latter scenario is more likely for two reasons. First, Africanized honeybee populations living in subtropical-temperate transitional zones of South America have retained a substantial component of European origin several decades after initial colonization (Lobo et al. 1989; Sheppard et al. 1991a; Lobo and Krieger 1992; Diniz et al. 2003; Quezada-Euán et al. 2003). Second, the WWR data of 2000 and 2001 indicated an influx of European alleles. In 2000, a strong heterozygote deficiency (P = 0.000) was accompanied by a reduction in the A. m. scutellata proportions, as shown by two of the three estimators depicted in Figure 4. In 2001, both estimators of nuclear admixture and the genetic distance measure (Figures 3 and 4) indicated continued reduction in the level of Africanization. In both years linkage disequilibrium was detected for the linked loci A28/ED1 (P < 0.05; data not shown). This trend was also mirrored by a reduction in the number of individuals assigned as “pure Brazilian” (Figure 8). It appears that when honeybees of largely African genetic composition began to compete for more limited resources, it was the more European component of the population that become favored. The reappearance of disequilibrium in 2000 and 2001 would be most parsimoniously explained by migration of honeybees with higher European proportions back into the WWR hybrid zone. Estimates of the genetic makeup of the new colonies appearing at WWR lend support to this assumption. In fact it suggests that movement into the zone likely occurred by both more strongly Africanized honeybees and honeybees with a more European component. This can be readily seen in Figure 8 in comparison with Figure 5.
Genetic composition of the post-Africanization population:
Within 5 years of detection of the first migrant colony of the A. m. scutellata mitotype, a panmitic European population was replaced by panmitic mixtures of African and European subspecies. As for feral neotropical populations, the WWR Africanized gene pool was predominantly of A. m. scutellata origin (mR = 62.8%; mY = 74.7%; mtDNA = 69%). However, both mitochondrial (Pinto et al. 2004) and nuclear European markers were more frequent in the WWR population than in most feral neotropical populations (Hall and Muralidharan 1989; Smith et al. 1989; Hall 1990; Sheppard et al. 1991b; Quezada-Euán and Hinsull 1995; McMichael and Hall 1996; Clarke et al. 2001, 2002; Hall and McMichael 2001).
The WWR post-Africanization gene pool was composed of a diverse array of recombinant classes with a substantial European genetic contribution. On the basis of data from the last two years of the study we expect that European genes will be maintained at moderate frequencies. As Hardy-Weinberg equilibrium has been largely maintained, the WWR population is best viewed as a “hybrid swarm” (Rinderer 1986; Rinderer et al. 1991). Given the current data, we suggest, as in Pinto et al. (2004), that this hybrid swarm represents what will likely become part of the more southern limit of the developing hybrid zone.
We are indebted to D. De Jong for providing the sample collection from Brazil; to A. Cavazos, R. Medrano, J. Teer, L. Drawe, T. Blankenship, and S. Glasscock for their help in the field work at the Welder Wildlife Refuge; to L. Ross for technical assistance in molecular biology; to J. B. Patton for writing the program used to simulate the hybrid populations; to E. Saillant and C. McKenna for assistance in the data analysis; to J. W. Bickham for graciously providing laboratory space for the genetic analysis. Microsatellite analyses were performed in the Center for Biosystematics and Biodiversity at Texas A&M University. The United States Department of Agriculture-Agriculture Research Service, Honey Bee Laboratory at Weslaco, Texas, and the Welder Wildlife Foundation were instrumental in the initialization and long-term support of many aspects of this research project. M. A. Pinto thanks Escola Superior Agrária de Bragança and the European Union program PRODEP II (Programa de Desenvolvimento Educativo para Portugal; Medida 5/Acção 5.3) for financial support. This is Welder Wildlife Foundation contribution no. 597. The Texas Legislative Initiative: Protection and Management of Honey Bees—Pollinators of Agricultural Crops, Orchards and Natural Landscapes funded this research.
Communicating editor: D. Begun
- Received August 16, 2004.
- Accepted April 14, 2005.
- Genetics Society of America