The plight of honey bees and other pollinators has gripped the attention of scientists and the general public alike, bringing increased scrutiny of the many environmental stressors on bees, particularly in managed agricultural habitats (1–4). It is widely held that the three main challenges pollinators face are parasites and pathogens (5–7), pesticides and other environmental contaminants (8–10), and nutritional deficiencies particularly as a result of habitat loss (6, 11–13). In addition, these challenges interact in critical and highly complex ways. This conventional wisdom, however, inherently rests upon the underlying genetic architecture of pollinator populations, which is rarely accounted for in studies of bee health. Since phenotypes are a function of genetics, the environment, and their interaction, quantifying the genetic background and relatedness of bees is imperative to fully understanding how environmental stressors may affect them.

In managed honey bees, population genetic structure has been investigated at several levels. On the evolutionary timeframe, phylogenetic studies have determined aspects of speciation and adaptive radiation (14–16). This provides an important historical perspective since it demonstrates the genetic diversity at the species or subspecies level, but since most imported Apis subspecies were of the M and C lineages (17), diversity at the evolutionary level does not likely impact how managed colonies interact with their local environment. At a more proximate and ecological timescale, Harpur et al. (18) tested managed honey bee populations in the US and Europe, then compared them to wild populations in Africa, Eastern Europe, and Western Europe. They concluded that there is significantly more admixture and genetic diversity in the introduced range, which prompted a debate (19, 20) about the merits and disadvantages of panmixis in managed honey bee populations (17, 21, 22). There have also been several studies comparing managed and putatively feral US honey bee colonies (23–26), showing differences in their relative genetic diversity. Finally, determining the genetic diversity of breeding programs has been of particular interest to researchers (27–29) due to a relatively small number of queen producers (mainly in Hawaii, California, and the southeast) accounting for the majority of genetic stock among beekeepers in the US (30).

Genetic diversity at these higher levels is critical to that at the colony level. One of the leading evolutionary explanations for hyperpolyandry (extreme female multiple mating) in honey bees is that the resulting increased intracolony genetic diversity is adaptively favored in many ways (reviewed in (31–34)). As can be seen in examples of genetic bottlenecking in other agricultural or livestock systems (35–40), it is important to the long-term sustainability of a healthy honey bee population to maximize genetic diversity and to avoid genetic bottlenecks among generations over the long term (30, 41–43).

Because of our previous work analyzing the genetic relatedness of individuals within colonies in an effort to estimate the effective paternity frequencies of honey bee queens (see below), we have a unique opportunity to conduct a meta-analysis on the individual genotypes of queens and their mates to quantify the degrees to which honey bees are related at the colony, operation, and population levels. Our goals were to quantify the average relatedness among breeding individuals in different honey bee populations, compare the pedigree relationships among breeding individuals, and determine the genetic structure of the US managed honey bee population.

We found significant differences in the average relatedness among the five breeding populations for all queen- and drone comparisons ( Figure 2 ). As expected, the average relatedness of queens to their mates (queen-mate) was not significantly different from zero for those in the AHB population (-0.035 ± 0.035), but those of all managed populations were significantly higher and not significantly different from each other (F _{4,207 =} 7.06, p<0.0001). Similarly, the AHB population had the lowest mate-mate relatedness (albeit significantly higher than zero; 0.049 ± 0.015) with all other populations being significantly higher but again not statistically different from each other (F _{4,207 =} 8.94, p<0.0001). The relatedness among queens within each population (queen-queen) followed a different pattern, with values in the two “control” populations as expected; queen relatedness in the AHB population were not significantly different from zero (-0.040 ± 0.033), and queens in the Experimental population had the highest relatedness (since they were known sisters; 0.361 ± 0.028). The three commercial populations were numerically intermediate and all above zero (F _{4,207 =} 64.4, p<0.0001), with the Breeding 1 population being statistically equivalent to the Experimental population, the Migratory population being significantly lower, and the Breeding 2 population being between the Migratory and AHB populations ( Figure 2 ).

Differences in average relatedness among drones mating with different queens within the same population demonstrated some particularly intriguing patterns ( Figure 2 ). Again, the drones in the AHB population were maximally outbred and had an average relatedness below zero (-0.032 ± 0.006). All other populations had average relatedness values above zero, with the Migratory and Breeding 2 population being statistically equivalent and significantly lower than the Breeding 1 population (F _{4,207 =} 51.9, p<0.0001). Interestingly, the Experimental population drone-drone relatedness was statistically lower than that of all commercial populations (0.016 ± 0.008), albeit still higher than the AHB population, suggesting that our experimental research apiary has a higher genetic diversity of drones than do the tested commercial populations.

Alternatively, when we estimated the key Mendelian relationships among the reproductives within and among colonies across the five populations ( Figure 3 ), we only found significant differences in the queen-queen comparisons (χ² ₁₂ ⁼ 94.6, p<0.0001) and none in the queen-mate (χ² ₈ ⁼ 0.01, p=1.00), mate-mate (χ² ₈ ⁼ 0.001, p=1.00), or drone-drone comparisons (χ² ₈ ⁼ 7.49, p<0.48). For queens, our AHB and Experimental populations were minimally and maximally inbred, respectively, as expected. The three commercial populations, however, were all intermediate, each with >48% of the queens being the genetic equivalent of supersisters (G=0.75), half-sisters (G=0.25), or cousins (G=0.125) ( Figure 3 ).

There were no differences in allelic richness among the six microsatellite loci ( Table 1 ), with all loci having 3.63-7.36 alleles in each study population (F _{5,24 =} 1.80, p=0.15). There were, however, significant differences in allelic richness across the studies (F _{4,25 =} 4.35, p<0.001), with the AHB population having significantly more alleles than the Breeding 1, Migratory, and Experimental groups with the Breeding 2 population being intermediate. The multilocus F_st values for all pairwise comparisons were highly significant ( Table 2 ).

The five experimental groups formed distinct clusters when analyzed with STRUCTURE ( Figure 4 ). When all data were analyzed together, the optimal number of distinct genetic populations was K = 4 based on (Ln P(D)) (54) and (Δ K) (55). The maximal value of Ln P(D) = -3718.5 for an L(K) of 6 based on six microsatellites from 211 queens. Groups 3 and 4 derive from the distinct AHB and Experimental populations, respectively, suggesting that there are two breeding populations in the three studies of managed stock ( Figure 4 ).

The results of our meta-analysis suggest that commercial breeding populations of managed honey bees are more inbred than maximally outcrossed but that there is sufficient and significant genetic variation depending on the level of analysis. While the five experimental datasets do not represent unique or separate breeding populations (i.e., some were surveys of multiple breeders whereas others were specific to certain operations) and thus do not represent a systematic sampling of all potential breeding stock in the US, they still represent many of the relatively small number of queen-producing operations and thus most managed honey bees. Since a small number of queen producers provide the majority of queens purchased by beekeepers each year (27, 30), our comparisons are valid representative examples of all potential breeding populations.

Our relatedness and pedigree analyses suggest that queens and drones within breeding populations can often be related by descent, thus demonstrating the potential for inbreeding. However, it appears that the overall effect on the average genetic relatedness within colonies is minimal. Perhaps most important is when queens mate with brothers or those closely related by descent, since this has the highest potential for resulting in homozygosity at the sex locus and other deleterious genetic combinations. We found that queen-mate relatedness among the five studies was significantly lower in the AHB feral population (as expected) compared to all others, suggesting that managed honey bees are not as outcrossed as they could be. However, all five groups had nearly identical queen-mate pedigrees, where ~70% of the queens mated with drones that were statistically unrelated (and less than 5% were genetically indistinct from brothers) irrespective of their mating population, including the unmanaged Africanized colonies. Another comparison of interest is the queen-queen analysis, which determines the likelihood of two queens being related by descent. This shows clear and intentional differences between the minimally inbred AHB and maximally inbred Experimental group. The three commercial groups were intermediate but much closer to having a higher likelihood of supersisters, half-sisters, or cousins. This finding is not surprising, given the standard methods of commercial queen rearing (56); while a given operation may use multiple breeder queens from which to graft young larvae to produce daughter queens, those from any given cohort likely derive from the same breeder queen and thus have higher average relatedness.

When looking at the population genetics of all queens, our results suggest that there is minimal structure within commercial honey bees, highlighting that the population in the US is not completely panmictic but nearly so (18). While accounting for the known AHB and Experimental populations, there appears to be only two distinct genetic populations across all queens from commercial managed stock. While our data do not permit any geographic inferences to be drawn, it is possible that this represents some genetic distinctions between queens reared in the southeast and west, where many if not most large-scale commercial queen producers are clustered. It would be interesting to conduct some further and systematic sampling and analysis of these two major queen-producing regions to quantify the degree of genetic overlap. Moreover, our results suggest that there is fairly clear, if limited, genetic introgression of AHB alleles into commercial stock (57–60).

This meta-analysis demonstrates that investigations into genetic diversity in honey bees can benefit from multiple perspectives. While population genetics can reveal overall trends and quantify the degree of potential diversity, Mendelian relationships are more critical at capturing functional bottlenecks that manifest at the colony level. Priority should be given to assessing genetic diversity within individual breeding populations at each of these levels to minimize inbreeding and hence optimize sustainability of the breeding stock of managed honey bees. In doing so, breeding programs should regularly test for allelic diversity (especially at the csd locus) to avoid inbreeding depression. Moreover, breeding programs with the goal of local environmental adaptation should eliminate or minimize admixture of genetic stock by maintaining large, closed breeding populations so that selection can manifest.

The data analyzed in this study is subject to the following licenses/restrictions: Previously published material. Requests to access these datasets should be directed to David Tarpy, drtarpy@ncsu.edu.

All authors equally contributed to the experimental design and final writing of the manuscript. JC was responsible for developing the RELATEDNESS simulations and pedigree comparisons, as well as sample processing and genotype generation in the original studies. DD was involved in sample collection, data processing, and population-genetic analyses, with DT supporting these efforts as well as writing the initial draft, statistical analyses, and figure creation. DT and DD were responsible for acquiring funding. All authors contributed to the article and approved the submitted version.