- Open Access
Methods to estimate breeding values in honey bees
© Brascamp and Bijma; licensee BioMed Central Ltd. 2014
- Received: 9 October 2013
- Accepted: 18 July 2014
- Published: 19 September 2014
Efficient methodologies based on animal models are widely used to estimate breeding values in farm animals. These methods are not applicable in honey bees because of their mode of reproduction. Observations are recorded on colonies, which consist of a single queen and thousands of workers that descended from the queen mated to 10 to 20 drones. Drones are haploid and sperms are copies of a drone’s genotype. As a consequence, Mendelian sampling terms of full-sibs are correlated, such that the covariance matrix of Mendelian sampling terms is not diagonal.
In this paper, we show how the numerator relationship matrix and its inverse can be obtained for honey bee populations. We present algorithms to derive the covariance matrix of Mendelian sampling terms that accounts for correlated terms. The resulting matrix is a block-diagonal matrix, with a small block for each full-sib family, and is easy to invert numerically. The method allows incorporating the within-colony distribution of progeny from drone-producing queens and drones, such that estimates of breeding values weigh information from relatives appropriately. Simulation shows that the resulting estimated breeding values are unbiased predictors of true breeding values. Benefits for response to selection, compared to an existing approximate method, appear to be limited (~5%). Benefits may however be greater when estimating genetic parameters.
This work shows how the relationship matrix and its inverse can be developed for honey bee populations, and used to estimate breeding values and variance components.
- Estimate Breeding Value
- Work Effect
- Virgin Queen
- Nuptial Flight
- Single Queen
Currently, honey bees (Apis mellifera) draw a lot of public and scientific attention because of increased colony losses ,, which are partly caused by infection with Varroa mites . Although selection is a promising way to improve Varroa tolerance of honey bees, estimation of breeding values is not common practice in this species ,. One reason is that it requires an organised collection of data on a relevant scale, which is rarely the case in honey bees. Currently, estimation of breeding values in honey bees is performed only in the German Beebreed program (http://www.beebreed.eu), for which breeding values are estimated from data that are collected annually on about 6000 colonies . For specifics on the genetic evaluation method used in the Beebreed program that we refer to as BER for Bienefeld, Ehrhardt and Reinhardt, please see .
Methodology for breeding value estimation in honey bees has drawn the attention of animal breeders -. They discussed the calculation of additive genetic relationships that account for the fact that the workers in a colony descend from a single diploid queen and 10 to 20 haploid drones. One approach that focused on the haplo-diploid nature of honey bees , suggested that an allelic relationship matrix that contains relationships between gametes instead of between individuals, can be adapted to the specifics of honey bee ancestry. Another approach focused on the uncertainty about the father of an individual  and suggested that methods developed for the use of mixed semen of sires can be adapted to honey bees. To our knowledge, these approaches have not been developed for implementation.
Breeding value estimation with an animal model builds on the work of Henderson , who derived the required inverse of the numerator relationship matrix using a decomposition of breeding values into Mendelian sampling terms. Because Mendelian sampling terms are mutually independent, the covariance matrix of these terms is diagonal, which facilitates inversion. However, it is not a diagonal matrix in honey bees , because the paternal contribution to the additive genetic relationship differs between workers in the same colony and workers in different colonies. Bienefeld et al.  solved this problem by using an approximation, in which both contributions are averaged and breeding values are estimated with an animal model. As a result, the matrix of Mendelian sampling terms is diagonal again, but the weighting of information of relatives is approximate.
The purpose of this paper was to develop a method, referred to as BB (for Brascamp and Bijma), to calculate the relationship matrix and its inverse for honey bee populations, in order to estimate breeding values and genetic parameters with an animal model. We used the approach of Henderson  as a starting point to derive the required procedures, taking in account the biology of the reproduction in honey bees. We also summarize the BER method and provide insight into the quantitative differences between the BB and BER methods, using Monte Carlo simulation in a simple example.
Reproduction of honey bees and colony observations
There are three types of individuals in honey bees: queens, workers and drones. Queens and workers are diploid, while drones are haploid. A colony of honey bees consists of a single fertilized queen, around ten thousand workers and several hundred drones. Workers contribute, for example, to the collection of pollen and nectar, the production of wax and nursing of the queen, but have no role in reproduction. Drones, in contrast, only serve for reproduction.
The description of the reproduction cycle in honey bees starts with a virgin queen. Soon after emerging from the brood cell, the virgin queen leaves the colony (nuptial flight) to mate in flight with multiple drones that come essentially from other colonies. These drones concentrate in so-called drone congregation areas, bringing together queens and drones in a range as large as 10 km. Drones die immediately after mating, which means that they can mate to a single queen only. Queens are mated only during their nuptial flight, or perhaps during a few nuptial flights within a small time slot and they cannot be mated again later in life. The queen stores a life lasting stock of millions of sperm cells in her spermatheca. After returning to their colony, mated queens produce two types of eggs, fertilized and unfertilized eggs. Fertilized eggs usually develop into diploid workers, while unfertilized eggs develop into haploid drones. Occasionally, an offspring of a fertilized egg receives a special diet from the workers and as a consequence develops into a virgin queen, which means that both workers and a virgin queen develop from a fertilized egg. The haploid drones that develop from unfertilized eggs have no father. They can be considered as flying gametes, and produce cloned sperm (i.e., all gametes produced by a drone are genetically identical).
Controlled mating of queens requires control of drones, which is possible only by restricting the presence of drone-producing queens with a particular pedigree on isolated mating stations (e.g. islands), or by artificial insemination. Under normal circumstances, in a colony, drones are produced along with workers, but the production of drones can be stimulated by management measures. Note that queens are always mated to multiple drones, both with natural mating and artificial insemination. Thus, the worker progeny of a queen descend from multiple drones. This situation resembles that with mixed semen in the case of e.g. pigs, for which the progeny of a sow derive from multiple boars. With respect to genetic relationships, the key difference between bees and mixed semen in pigs is that each piglet descends from a genetically unique paternal gamete, while subsets of the workers in a colony descend from the same drone and therefore from genetically identical paternal gametes.
Drones cannot be not traced and it is unknown how many and which drones have mated to the queen. As a consequence, the contribution of each drone-producing queen to the offspring of the queen is unknown. For this reason, the group of drone-producing queens can be treated as a single “individual”, which we will refer to as the sire of the workers of the colony. In Figure 1, for example, the three drone-producing queens in 1b together constitute the sire of the workers in the colony of queen 1a. By grouping the drone-producing queens into a single sire, each individual in the pedigree has precisely two parents, a queen and a sire. This grouping makes it easier to trace the pedigree without loss of information.
Observations in honey bees are on colony performance, and may relate to traits like honey production, behaviour and disease resistance . The performance of a colony is affected by the joint genetic effects of the ten thousand workers (called worker effect) and by the genetic effects of the queen (called queen effect). Colony performance results from the action of the workers and the interaction between workers, but also from the effects of the queen on the workers, for example, due to the number of workers produced or by producing pheromones that affect worker behaviour. However, workers also affect the behaviour of the queen. Despite these different interactions, the performance of a colony can be partitioned into an additive worker effect and an additive queen effect, based on the principle of least squares . Conceptually, this is similar to defining the average effect of an allele for a locus showing dominance, and to maternal effects in mammals . Several studies , have shown that the contribution of queen effects to colony performance is considerable, although smaller than that of worker effects, while the genetic correlation between worker and queen effects is negative.
In the following paragraphs, we consider three types of individuals: (i) queens, (ii) sires, and (iii) groups of workers in a colony, referred to as worker groups. Queens are single individuals, while sires and worker groups are aggregates of individuals. With this categorisation, we cover individuals responsible for the phenotypes (queens and worker groups) and individuals in the pedigree (queens and sires). To emphasize that worker groups and sires consist of groups of individuals, we will write their breeding values as averages, using , while using A for the breeding value of a queen. Since breeding values are to be estimated on all three categories, the size of the numerator relationship matrix will be twice the number of queens (because each colony has one queen and one worker group) plus the number of sires.
Here and are the additive genetic variances for worker and queen effect, respectively, and r G is the genetic correlation between these effects.
In the next section, we develop the method to derive A− 1 that is needed in Equation (4).
Numerator relationship matrix
Equations (6) through (8) hold for honey bees as well, but D is no longer a diagonal matrix. In the following, we derive the diagonals and off-diagonals of D, considering the three types of individuals defined above: queens, sires and worker groups. Because D is the same for all traits of interest, we do not distinguish between worker and queen effects, and therefore drop the W and Q superscripts.
Diagonal elements of D
The first component in Equation (13) represents the variance due to the Mendelian sampling of maternal gametes, the second one the variance due to the Mendelian sampling of gametes of an individual drone-producing queen, and the third one the variance due to the sampling among drone-producing queens. Note that this equation can be applied when both parents of individual i are known. If this is not the case, refer to Appendix 1.
Hence, Equation (19) shows that the worker group has a non-zero sampling term merely because workers may descend from the same drone; otherwise would average to zero. Finally, diagonal elements for worker groups follow from Equation (18).
Off-diagonal elements of D
Covariance between sampling terms of full-sibs cov(δ FS )
In honey bees, full-sibs are the offspring of the mating between a queen and a sire. Within a colony, some pairs of workers are full-sibs in the ordinary sense (when they descend from a common queen and a common drone-producing queen, but from different drones) with an additive genetic relationship of , ignoring inbreeding. A pair may also descend from the same drone, resulting in , or from different drone-producing queens, resulting in .
Usually in animal breeding, the covariance between breeding values of relatives is fully accounted for by the pedigree. In general, however, this requires two conditions. The first condition is that, conditional on the pedigree, Mendelian sampling terms of offspring are independent. In honey bees this is not the case for full-sibs, because they may descend from the same drone, in which case their paternal Mendelian sampling terms are identical. The second condition is that the pedigree fully accounts for the contributions of parents to offspring. Usually in animal breeding, this condition is met, because a parent contributes precisely half the genes of an offspring. However in a honey bee pedigree, this condition is not met because the sire is an aggregate of multiple drone-producing queens and the contribution of individual drone-producing queens to offspring varies among the drone-producing queens that constitute a sire. This will occur by chance, even when the a priori expected contribution is the same for all drone-producing queens that make-up a sire but the pedigree accounts for only the average contribution of a drone-producing queen to the offspring, which is given by the in Equations (9) and (14). Variation among drone-producing queens in their contribution to offspring creates a paternal covariance among full-sibs that exceeds the that is accounted for by the pedigree, and thus creates a covariance between the δ terms of full-sibs. Thus, the δ terms of full-sibs may be correlated because (i) sibs may descend from the same drone, and (ii) the contribution to offspring may vary among drone-producing queens.
The first term in equation (21) arises from the probability that two full-sibs descend from the same drone. The second term arises from variation in the contribution of individual drone-producing queens to offspring, around the expected value of . Thus, in the second term, the term represents subtraction of the covariance already accounted for by the pedigree.
Finally, off-diagonals of D are obtained from substituting Equation (21a) into Equation (20). Thus, when the number of offspring per parent follows a Poisson distribution, the covariance between Mendelian sampling terms of full-sibs depends only on the relatedness between drone-producing queens (a ss ) and on the number of drones mated to a queen. Under the assumption that the number of drones of a single drone-producing queen that mates to the queen follows a Poisson distribution, there is no covariance between sampling deviations of paternal half sibs (i.e., between two offspring of the same sire but of a different queen; see Discussion). Whether the assumption of a Poisson distribution is realistic will be addressed in the Discussion.
Note that the BER method does not implement off-diagonal elements in D (see below); here, we merely present Equation (21b) to show the outcome of cov(δ FS ) for the p2 proposed by . Note that for S-1 approaching S, Equation (21b) approaches (21a).
Construction of D and D −1
In this equation, the first term represents the additive genetic relationship between drone-producing queens because they descend from the same dam, the second term relates to the case when they descend from the same drone, which has probability p1 and a paternal relatedness of , the third term relates to the case when they descend from the same drone-producing queen but from a different drone, which has probability (p2 − p1) and a paternal relatedness of , the fourth term relates to the case when they descend from different drone-producing queens, which has probability (1 − p2) and a paternal relatedness of , and the last term accounts for the additive genetic relationship between dam and sire of the drone-producing queens.
As a result, D is a block-diagonal matrix, each block representing the offspring of a single queen, i.e. the combination of a queen and a sire. Chronologically, such a block starts with a single individual, being the worker group that descends from that queen. When the queen is selected to breed new queens, the queens in its progeny will be added to the block. Moreover, when the queen is selected to breed drone-producing queens, then one or more sires will be added to the block. The size of a block, therefore, equals 1 plus the number of queens plus the number of sires that descend from the mother queen. Thus, a block contains at maximum three distinct diagonal values, one for the worker group, one for queens, and one for sires. All off-diagonals within a block are equal, and equal to the diagonal element for the worker group (Equation (20)). Off-diagonals outside blocks are 0.
Since D is a block-diagonal matrix, the inverse of D is also a block-diagonal matrix, each block being the inverse of the corresponding block of D. Since blocks of D have a specific structure, with at maximum three distinct values, D− 1 can be obtained analytically, e.g., with the help of equation-solving software such as Mathematica . However, since blocks of D can have different numbers of queens and sires, there are multiple analytical solutions, each of which is a complicated expression. Therefore, since the size of the blocks is usually small, numerical inversion of each block is easy and more practical and, thus, we do not present the analytical inversion of D here.
The Bienefeld, Ehrhardt and Reinhardt (BER) method 5
In the BER method, both the sire and the worker group are treated as single individuals, not as virtual individuals that consist of a group of individuals. As a consequence, the diagonal elements of D are neither affected by the number of drone-producing queens nor by the number of drones involved in mating, so that the sampling variances are computed in the same manner as those of queens.
The purpose of the simulation was to study properties of estimated breeding values from the BER and BB methods, and to compare the estimated breeding values. In the simulation, breeding values are generated for three types of individuals, queens, sires and worker groups, and subsequently the phenotypes for colonies are generated.
Simplified selection cycle in the honey bee selection programme
Selection of queens to produce drone-producing queens (sires); birth of sires
t + 1
Selection of queens to produce full-sib groups of queens to be mated to sires
Use of sires for mating to groups of queens, each group being the progeny of a selected queen
t + 2
Test of colonies of the set of full-sibs born in t + 1
t + 3
Selection of queens and birth of next-generation queens to be mated to sires born in t + 2
Selection of queens to produce drone-producing queens
The basic simulation starts with NSY (number of sires per year) base sires generated in years 1, 2 and 3, and NQY (number of queens per year) base queens generated in years 2 and 3. The base queens were simulated as full-sib groups with NQ (number of full-sibs per group) individuals each because that is also the structure in future generations. Sires from years 1 and 2 are mated randomly to queens in years 2 and 3 with an equal number of NFQ mates for each sire, each producing NQ full-sibs. From year 4, the NSY queens with the highest estimated breeding value according to Equation (2) are selected to produce the next generation of sires. The NSYxNFQ queens with the highest estimated breeding value are selected to produce NQ queens each. Allocation of mates (sires) to queens is random, but each queen within a full-sib group is mated to the same sire and each sire is mated to the same number of queens. The queens that will produce sires are therefore selected NSY out of NQY, and the queens that will produce queens are selected 1 out of NQ.
Two breeding values are simulated for each individual, a breeding value for the worker effect and a breeding value for the queen effect. To allow for a correlation between these two breeding values, random samples for an individual are taken from a bivariate normal distribution (using the function mvrnorm from the package MASS in R; http://cran.r-project.org/package=MASS). Because Mendelian sampling terms are correlated between offspring of the same queen, sampling terms were constructed as the sum of two independent components: a component specific to each individual and a component common to all offspring of the same queen.
Using the simulated data and pedigree, the inverse numerator relationship matrix A− 1 (Equation (8)) was created using either method BB or method BER, and breeding values were estimated by solving Equation (4). The criterion to select queens to produce the next generation of queens and sires for method BB is given by Equation (2). For method BER, the factor of for the estimated breeding value of the sire in Equation (2) is replaced by q.
Parameter values used in the simulation were = 1, = 0.5, = 2 and r G = −0.5, which are in line with estimates reported by . Based on , the number of drone-producing queens that constitute a sire (S) was equal to 8 and the number of drones mating to a queen (D) was equal to 12. We analysed the simulated data for a small example, in which the only fixed effect was the mean. NSY was equal to 5, the number of full-sib groups to which a sire is mated (NFQ) was equal to 5 and NQ was fixed at 3. We simulated 20 years of data, including colony performance of queens born in year 20.
First, the properties of estimated breeding values (EBV) were investigated using 1000 replicated schemes without selection. The quality of EBV was judged by the regression coefficient of the true (i.e., simulated) breeding value (TBV) on the EBV and by the correlation coefficient between TBV and EBV in year 20. We chose those criteria because the regression coefficient of TBV on EBV should be equal to 1 with BLUP (best linear unbiased prediction), while response to selection on EBV is proportional to the correlation coefficient. We did not implement the correction factor from Equation (30), since this does not affect results because regression and correlation coefficients were calculated within one generation. Second, we compared response to selection between the two methods, again using 1000 replicates. Because selection took place within years (non-overlapping generations), we did not implement the correction factor of Equation (30) here either, since it did not affect results.
Regression coefficients of true breeding values on estimated breeding values 1 obtained from two methods (BB and BER)
Correlations between the true and estimated breeding values and cumulative responses to selection with two breeding value estimation methods (BB and BER) 1
Cumulative response to selection3
Responses to selection are also in Table 3. The annual responses to selection started slowly in early years and remained somewhat irregular in later years. There was a strong similarity in results for the two methods in that respect. This irregularity is caused by the structure of the simulation. The simulation started with the simulation of base sires in years 1 to 3 and base queens in years 2 and 3. A first batch of offspring (born in year 4) is produced from base sires of year 1 and base queens of year 2 and a second batch (born in year 5) from base sires of year 2 and base queens of year 3. Genetically, progeny of these batches of offspring mixed in later years due to the fact that the dams of the sires were three years old at birth of the next generation of queens, while the dams were two years old at birth of the next generation of queens. However, this mixing developed fairly slowly and was delayed by the fact that pairs of queens and sires were selected to produce the next generation. Similar results have been observed in simulations of breeding schemes with overlapping generations in dairy cattle . The cumulative selection response differed little between the two methods BB and BER but were 5% higher with method BB compared to method BER in years 8 to 10, and 4% higher in years 18 to 20.
In this paper, we derived a method to calculate the relationship matrix and its inverse for honey bee populations, which is required to estimate breeding values and genetic parameters. The situation in honey bees differs from the usual situation in farm animal breeding, because of the honey bees’ mode of reproduction. The first major difference is that two full-sibs may carry identical paternal gametes. This occurs because sires (drone-producing queens) produce drones which may be considered as flying gametes that produce many identical sperm cells. Because a drone can mate to a single queen only, paternal half-sibs always carry different paternal gametes. Consequently, the paternal contribution to the additive genetic relationship between full-sibs differs from that between half-sibs, which results in a block diagonal D matrix of covariances between Mendelian sampling terms. Off-diagonals of those blocks equal the covariance between sampling terms of full-sibs. The second difference is that selection candidates (queens) are mated early in life, before they can be selected as parents. As a consequence, selection is not of individual dams but of matings from which breeding stock can be produced after the estimation of breeding values. Thus, the selection target is the breeding value of a future queen from this mating, which equals half the breeding values of both mates plus the part of Mendelian sampling that is common to all progeny of these mates. This also implies that the EBV of such a future queen equals the EBV of the colony. Another difference from traditional animal breeding is that the “father” of a queen is usually unknown because the drones that mate with the queen come from multiple drone-producing queens. In this context, our work follows that of Dempfle , who discussed the consequences of mixed semen for the estimation of breeding values, rather than focussing on the haploid nature of the drones ,.
Equation (21) gives a general expression for the covariance of the Mendelian sampling terms of full-sibs. This covariance depends on the variance of the number of offspring per drone (Equation (22)) and the variance of the number of drones per drone-producing queen (Equation (23)). Without further knowledge, a Poisson distribution of family sizes is a common choice, which leads to Equation (21a). Numerically, this equation differs very little from Equation (21b), which results from substituting the probabilities from the BER method (Equation (23b)) into Equation (21), but these probabilities do not have a theoretical basis. However in reality, the assumption of a Poisson distribution of family sizes does not seem to hold, since a review of the literature  suggests that the proportion of progeny descending from different drones deviates from Poisson. Furthermore, results of an experiment using drone-producing colonies each producing a similar number of drones, suggested that drone-producing queens that contribute a higher proportion of drones to matings also produce drones with a higher proportion of offspring in colonies . The Poisson distribution arises when variation in contributions is entirely by chance, i.e., when a priori the expected number of offspring is equal for each drone-producing queen and for each drone. When there are systematic differences between drone-producing queens or drones in the expected number of offspring, then the variance in contributions will be larger than in the case of a Poisson distribution, which implies a larger covariance between sampling terms of full-sibs. Numerically, this effect may be neutralised by assuming a smaller number of drones, i.e., by using an effective number of drones rather than the actual number of drones. Note that in this context,  used a number of drones equal to 12, while a consensus number is around 16 . When assuming a Poisson distribution, the covariance between Mendelian sampling terms for half sibs is 0. This is, however, not true if some drone-producing queens are systematically more successful to contribute drones to matings .
In practice, in the Beebreed program, inbreeding coefficients are computed for possible planned matings that are not yet included in the pedigree (http:www.beebreed.eu). Efficient methods to compute inbreeding coefficients have been derived ,, based on . These methods exploit the fact that the D matrix is a diagonal matrix. We derived a modification to , which takes into account the fact that the D matrix contains off-diagonal elements in the bee breeding case. This method, however, requires the whole pedigree of an individual to be searched for the occurrence of parents that are full-sibs, which may be very time-consuming. As an alternative, the A matrix of the pedigree may be kept in memory such that the required inbreeding coefficients can be used in Equation (25).
In the development of the methods and analyses presented here, we used the current mating system applied in the Beebreed program as a starting point. This implies that drone-producing queens are full-sibs from a shared dam and sire. That may not be the case for parts of the pedigree or for other programs. For those cases, we suggest to include the individual drone-producing queens in the pedigree, with diagonal elements in D equal to those of individual queens, combining Equations (13) and (10). Elements in M then need to be adapted to reflect the fractions that are contributed by these individual drone-producing queens. Without prior knowledge on these fractions, we suggest to use equal fractions as an approximation, although this may not be true in reality ,.
We have presented methodology to construct the relationship matrix and its inverse for honey bee populations, which is required in the mixed model equations used for the estimation of breeding values and genetic parameters. The method allows for different assumptions on the contribution of drones and drone-producing queens to offspring, and is exact if those assumptions are correct. The method yields EBV that are unbiased predictors of TBV. We also carried out an exploratory comparison with the BER method  that is currently used in practice and weighs information on relatives, differently. Although EBV obtained with the BER method were biased, selection candidates were ranked similar to those of our method and the response to selection was only slightly lower than with our method. This suggests that suboptimal weighting of information from relatives has limited impact on the ranking of selection candidates. It remains to be seen whether this conclusion extends to the estimation of genetic parameters.
EWB had the initial idea of looking into the estimation of genetic parameters in honey bees based on the theoretical framework of Bienefeld et al. for breeding value estimation. PB suggested focusing on this theoretical framework which led to this paper with a focus on breeding value estimation. EWB and PB jointly developed the theory and wrote the paper, with EWB writing the first draft. EWB designed the simulation and programmed it in R. Both authors read and approved the final manuscript.
Variance of Mendelian sampling terms
In this Appendix, first we derive Equation (13), the Mendelian sampling variance for queens with known parents for method BB. Subsequently, we derive the variances when parents are unknown, for queens and sires. Finally, we repeat all steps for method BER. Equation numbers between brackets refer to equations in the main text.
A queen with unknown parents (base queen). In that case, A i = δ i and so .
A sire with unknown parents (base sire). In that case, and , because the drone-producing queens constituting the sire are assumed to be unrelated.
An individual with a known dam but an unknown sire. In that case, and such that .
An individual with an unknown dam, but a known sire. Then, and and after some rearrangement: .
A sire with an unknown dam and a known sire. This is very unlikely because it implies that a sire is bred from a dam for which neither pedigree nor breeding values are available. Nevertheless, in that case and , as when both parents are known with δ i as under point 4 of this paragraph, such that where var(δ i ) is as under point 4 and as before, when both parents are known. Therefore, .
A sire with a known dam and an unknown sire. This also seems odd. Nevertheless, in that case , and with δ i as under point 3. Now, , with var(δ i ) as under point 3. Covariances between δ i and δ j are equal to 0 because no covariance due to drones is involved.
When neither the dam nor the sire is known .
In the case when only the dam is known, and . In the case when only the sire is known, A i = qA s + δ i , and .
Here, we derive the probability p1 that two offspring of a queen descend from the same drone, and the probability p2 that two offspring descend from the same drone-producing queen, under the assumption that the number of drones per drone-producing queen and the number of offspring per drone follow a Poisson distribution. Take the average number of drones per drone-producing queen to be equal to D S and the average number of offspring per drone to be equal to N D . Furthermore, take the total number of offspring per sire to be equal to T, such that T = SD S N D and the number of offspring per drone-producing queen to be equal to N S .
The proportion of offspring that derives from a particular drone is .
It follows that .
The average number of drones per drone-producing queen equals . This number follows a Poisson distribution such that .
We have defined c S as the proportion of offspring from a particular drone-producing queen, where on average , using T = SD S N D .
The contribution of JWM Bastiaansen to write an effective selection algorithm in R is gratefully acknowledged. PB acknowledges financial support of the Netherlands Organization for Scientific Research (STW-NWO).
- Ghazoul J: Buzziness as usual? Questioning the global pollination crisis. Trends Ecol Evol. 2005, 20: 367-373. 10.1016/j.tree.2005.04.026.View ArticlePubMedGoogle Scholar
- Van der Zee R, Pisa L, Andonov S, Brodschneider R, Charrière JD, Chlebo R, Coffey MF, Crailsheim K, Dahle B, Gajda A, Gray A, Drazic MM, Higes M, Kauko L, Kence A, Kence M, Kezic N, Kiprijanovska H, Kralj J, Kristiansen P, Martin Hernandez R, Mutinelli F, Nguyen BK, Otten C, Özkırım A, Pernal SF, Peterson M, Ramsay G, Santrac V, Soroker V: Managed honey bee colony losses in Canada, China, Europe, Israel and Turkey, for the winters of 2008–9 and 2009–10. J Apicult Res. 2012, 51: 100-114. 10.3896/IBRA.18.104.22.168.View ArticleGoogle Scholar
- Rosenkranz P, Aumeier P, Ziegelmann B: Biology and control of Varroa destructor. J Invertebr Pathol. 2010, 103: 96-119. 10.1016/j.jip.2009.07.016.View ArticleGoogle Scholar
- Büchler R, Andronov S, Bienefeld K, Costa C, Hatjina F, Kezic N, Kryger P, Spivak M, Uzunov A, Wilde J: Standard methods for rearing and selection for Apis mellifera queens. J Apicult Res. 2013, 52: 1-30. 10.3896/IBRA.1.52.1.07.View ArticleGoogle Scholar
- Bienefeld K, Ehrhardt K, Reinhardt F: Genetic evaluation in the honey bee considering queen and worker effects – A BLUP-animal model approach. Apidologie. 2007, 38: 77-85. 10.1051/apido:2006050.View ArticleGoogle Scholar
- Smith SP, Allaire FR: Efficient selection rules to increase non-linear merit: application in mate selection. Genet Sel Evol. 1985, 17: 387-406. 10.1186/1297-9686-17-3-387.PubMed CentralView ArticlePubMedGoogle Scholar
- Tier B, Sölkner J: Analysing gametic variation with an animal model. Theor Appl Genet. 1992, 85: 868-872.Google Scholar
- Dempfle L: Problems in the use of the relationship matrix in animal breeding. Advances in Statistical Methods for Genetic Improvement in Livestock. Edited by: Gianola D, Hammond K. 1990, Springer Verlag, Berlin, 454-473. 10.1007/978-3-642-74487-7_20.View ArticleGoogle Scholar
- Henderson CR: A simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values. Biometrics. 1976, 32: 69-83. 10.2307/2529339.View ArticleGoogle Scholar
- Bienefeld K, Pirchner F: Heritabilities for several colony traits in the honeybee (Apis mellifera carnica). Apidologie. 1990, 21: 175-183. 10.1051/apido:19900302.View ArticleGoogle Scholar
- Fisher RA: The correlation between relatives on the supposition of Mendelian inheritance. Philos Tran R Soc Edin. 1918, 52: 399-433. 10.1017/S0080456800012163.View ArticleGoogle Scholar
- Willham RL: The covariance between relatives for characters composed of components contributed by related individuals. Biometrics. 1976, 19: 18-27. 10.2307/2527570.View ArticleGoogle Scholar
- Ehrhardt K, Böchler R, Bienefeld K: Genetic parameters of new traits to improve the tolerance of honeybees to varroa mites. In Proceedings of the 9 th World Congress on Genetics Applied to Livestock Production: 1–6 August 2010; Leipzig. 2010 , [http://www.kongressband.de/wcgalp2010/assets/pdf/0565.pdf]
- Henderson CR: Sire evaluation and genetic trends. J Anim Sci. 1973, 1973: 10-41.Google Scholar
- Quaas RL: Computing the diagonal elements and inverse of a large numerator relationship matrix. Biometrics. 1976, 32: 949-953. 10.2307/2529279.View ArticleGoogle Scholar
- Wolfram Mathematica., , [http://www.wolfram.com/mathematica/]
- Ducrocq V, Quaas RL: Prediction of genetic response to truncation selection across generations. J Dairy Sci. 1988, 71: 2543-2553. 10.3168/jds.S0022-0302(88)79843-4.View ArticleGoogle Scholar
- Schlüns H, Moritz RFA, Lattorff MG, Koeniger G: Paternity skew in seven species of honeybees (Hymenoptera: Apidae: Apis). Apidologie. 2005, 36: 201-209. 10.1051/apido:2005006.View ArticleGoogle Scholar
- Kraus FB, Neumann P, Scharpenberg H, Van Praagh J, Moritz RFA: Male fitness of honeybee colonies (Apis mellifera L.). J Evol Biol. 2003, 16: 914-920. 10.1046/j.1420-9101.2003.00593.x.View ArticlePubMedGoogle Scholar
- Luo Z: Computing inbreeding coefficients in large populations. Genet Sel Evol. 1992, 24: 305-313. 10.1186/1297-9686-24-4-305.PubMed CentralView ArticleGoogle Scholar
- Sargolzaei M, Iwaisaki H, Colleau J-J: A fast algorithm for computing inbreeding coefficients in large populations. J Anim Breed Genet. 2005, 122: 325-331. 10.1111/j.1439-0388.2005.00538.x.View ArticlePubMedGoogle Scholar
- Matilla HR, Seely TD: Extreme polyandry improves a honey bee colony’s ability to track dynamic foraging opportunities via greater activity of inspecting bees. Apidologie. 2014, 45: 347-363. 10.1007/s13592-013-0252-3.View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.