Genetic diversity and population structure of two subspecies of western honey bees (Apis mellifera L.) in the Republic of South Africa as revealed by microsatellite genotyping

Honey bee sampling and DNA extraction

In April/May 2013 and April 2014, a total of 60+ bees were collected from each of over 1,100 managed honey bee colonies located in 75 different apiaries grouped into 29 geographical regions (1–3 apiaries/region) across the A.m. scutellata and A.m. capensis distribution in RSA (Mortensen et al., 2016). A region was defined as a cluster of 1–3 apiaries that were located close to one another (i.e., ∼ 5–50 km between apiaries) or was otherwise isolated by itself in instances where a region was composed of only one apiary. The regions were named for the closest city, town, village, or other geographically defined, nearby area. Detailed information on the sampling regions and apiaries, the number of bees analyzed per apiary and the geographical coordinates of the sampling regions are shown in Table 1.

The bees were collected and stored by colony in a 50 ml tube containing ≥95% ethanol. The samples were transported (per USDA-APHIS regulations) to the Honey Bee Research and Extension Laboratory at the University of Florida where the identity of each colony (A.m. capensis, A.m. scutellata or hybrid) was determined using classical morphometrics conducted on at least ten bees from each colony (for procedure, Bustamante, Baiser & Ellis, in press). The results of the morphometric classification are shown in Table S1. Colony identity was cross-checked using a distribution map published by Hepburn & Radloff (1998). For later analysis, the sample tubes of bees were kept at −80 °C.

For this study, 813 worker bees were sampled from 101 selected colonies across the 75 different apiaries, with 22–33 bees used per region (5–10 bees per colony). The DNA was extracted from the dissected thoraces of each individual bee using Wizard^® Genomic DNA Purification kit (Promega, Durham, NC, USA) according to the manufacturer’s instructions. DNA quality was assessed using a 1% agarose gel and quantified using Qubit^® 3.0 Fluorometer based on manufacturer’s guidelines (Thermo Scientific Inc., Waltham, MA, USA).

Microsatellite genotyping

Nineteen microsatellite loci were selected from previous studies and deemed appropriate candidates for the genotyping of the target bee populations (Shaibi, Lattorff & Moritz, 2008; Chen et al., 2005; Alburaki et al., 2013; Rowe et al., 1997; Solignac et al., 2003; Weinstock et al., 2006). The forward primers were labeled at the 5′ end with one of four fluorescent dyes (FAM, VIC, PET and NED; Applied Biosystems, Foster, CA, USA) and distributed into six multiplex PCR reactions. A list of labeled primers and associated information is included in Table S2. PCR reactions were performed in a thermal cycler (Eppendorf, Hamburg, Germany) in 20 µl volumes containing 10 µl of Master Mix Maxima Hot Start 2X Qiagen, 0.6 µl of each primer (10 pmol/µl) and one µl of DNA (>10 ng/µl). The PCR program started with a denaturation phase at 95 °C for 5 min, followed by 35 cycles of 30 s at 95 °C, 30 s at 55 °C, 30 s at 72 °C and a final elongation phase at 72 °C for 30 min. The resulting PCR products were diluted 25–100 fold before genotyping and submitted to the University of Florida’s Interdisciplinary Center for Biotechnology Research for sequencing. A volume of 990 µl Hi-Di™ Formamide (Catalog # 4311320; Life Technologies, Grand Island, NY, USA) was mixed with 10 µl of the GeneScan™ 600 LIZ size standard (Catalog # 1406056; Life Technologies, Grand Island, NY, USA). A volume of 10 µl of the solution was mixed with one µl of PCR product per sample and distributed in a 96 well plate. Fragment lengths were determined using a 3730*1 DNA Analyzer (Model 325-0020; Applied Biosystems, Foster City, CA, USA) and alleles were scored using GeneMarker version 2.4.0 (Applied Biosystems, Foster City, CA, USA) (Hulce et al., 2011). All electropherograms were checked manually and confirmed for further statistical analyses.

Statistical analyses

Genetic diversity

In an analysis of 19 microsatellite loci, we detected a high level of genetic variability based on the mean number of alleles observed per microsatellite locus and heterozygosity within and among examined honey bee regions. The high level of genetic diversity observed for all honey bee groups was consistent with that shown in previous studies (Estoup et al., 1995; Franck et al., 1998; Jaffe et al., 2010; Eimanifar et al., 2018a, 2018b). This may be true in RSA given the high densities of wild honey bee colonies in certain areas of the country (Vaudo et al., 2012a, 2012b).

The microsatellites we used tended to classify bees from certain regions as hybrids, more so than did SNPs, which tended to separate the subspecies more predictably (Eimanifar et al., 2018b). In fact, both the A.m. capensis and A.m. scutellata regions in RSA were smaller when defined by microsatellites than they were when defined using SNPs or whole mitogenomes (Eimanifar et al., 2018a, 2018b). Consequently, we feel the microsatellites possibly overestimated the size of the hybrid zone. In contrast, it is possible that the microsatellites provided evidence of greater hybridization between the two subspecies than first thought.

As expected, honey bees in the hybrid region in RSA contain a gene pool derived from hybridization of A.m. scutellata and A.m. capensis bees. The hybrid bees in our study had many polymorphic alleles that likely originated from recombination of alleles of the two subspecies (Alda & Doadrio, 2014). We suspect that the occurrence of new alleles in the hybrid bees could be a result of hybridization of the two subspecies, as is known to happen for populations in hybrid zones in general (Arnold et al., 1999).

The mean number of alleles observed across all loci (Na = 10.01) and the mean value of observed heterozygosity (H_obs = 0.75) were greater for bee populations in our study (N = 716) than those reported in the literature for Algeria (N = 414, Na = 9.5, H_obs = 0.69) (Loucif-Ayad et al., 2015), Croatia (N = 225, Na = 9.25, H_obs = 0.67) (Muñoz et al., 2009) and on Rodrigues Island (N = 524, Na = 7.63, H_obs = 0.64) (Techer et al., 2016). The greater mean number of alleles/locus in our study than in others likely resulted from our large sample size and generally large populations of the two subspecies we examined (Loucif-Ayad et al., 2015). African honey bees exhibit several exceptional characteristics, including pronounced migratory behavior (absconding and swarming), which possibly could explain their high level of genetic diversity with low levels of differentiation (Franck et al., 1998; Hepburn & Radloff, 1998).

Population genetic structure

The set of microsatellite markers we used yielded a distinctive population structure within the examined groups, though the microsatellite results did not align completely with those of the morphometric (Bustamante, Baiser & Ellis, in press) and other molecular results (Eimanifar et al., 2018b). Nevertheless, our STRUCTURE results support the presence of some genetic structure at the subspecies level, reinforcing other work showing that there are two different gene pools existing in A.m. capensis and A.m. scutellata honey bees from the RSA (Eimanifar et al., 2018a, 2018b).

We found additional evidence of the allelic introgression of A.m. capensis from their native range into that of A.m. scutellata regions (Neumann & Moritz, 2002; Moritz et al., 2008). These results highlight the decades-old, continuing impact of the “capensis calamity” on A.m. scutellata colonies in the RSA. Our STRUCUTRE results (Fig. 3) suggest that the A.m. capensis/scutellata hybrid zone has moved north and east from that originally delineated by Hepburn & Radloff (1998) over two decades ago.

In conclusion, honey bees in the RSA are composed of two genetically different populations, both with significant genetic structure. Given our large sample size and analytical tools, we show two distinct evolutionary units, though the results did not match those of earlier morphometric and molecular analyses (Eimanifar et al., 2018b; Bustamante, Baiser & Ellis, in press), suggesting that the microsatellites we tested were not sufficient for subspecies identification purposes, especially for Cape and hybrid bees. Nevertheless, our data highlight the considerable genetic diversity within both populations, a possibly larger-than-expected hybridization zone between the natural distributions of A.m. capensis and A.m. scutellata, and evidence of the expansion of A.m. capensis into A.m. scutellata regions. The genetic introgression of A.m. capensis into A.m. scutellata territory continues from the days of the “capensis calamity” and highlights the invasive potential of A.m. capensis. Our data suggest that the hybrid zone has expanded in the RSA beyond that shown using morphological data (Hepburn & Radloff, 1998; Bustamante, Baiser & Ellis, in press), possibly owing to migratory beekeeping activities in the RSA.