Honey yield

This section shows how to simulate one colony trait, honey yield, that is influenced by the queen and workers as well as the environment. We will achieve this by:

setting base population quantitative genetic parameters,
inspecting individual values in the base population,
inspecting individual values in a colony,
calculating colony value,
calculating multi-colony values, and
selecting on colony values.

Base population quantitative genetic parameters

AlphaSimR, and hence SIMplyBee, simulates each individual with its corresponding genome, and quantitative genetic and phenotypic values. To enable this simulation, we must set base population quantitative genetic parameters for the traits of interest in the global simulation parameters via SimParamBee. We must set:

the number of traits,
the number of quantitative trait loci (QTL) that affect the traits,
the distribution of QTL effects,
trait means, and
trait genetic and environmental variances - if we simulate multiple traits, we must also specify genetic and environmental covariances between the traits.

In honeybees, the majority of traits are influenced by the queen and workers. There are many biological mechanisms for these queen and workers effects. Depending on which caste is the main driver of the trait (the queen or workers), we also talk about direct and indirect effects. For example, for honey yield, workers directly affect honey yield by foraging, while the queen indirectly affects honey yield by stimulating workers via pheromone production. The queen and workers effects for a trait can be genetically and environmentally independent or correlated (usually negatively).

Here, we will simulate two traits to represent the queen and workers effects on honey yield. From this point onward we will use the terms the queen effect and queen trait interchangeably. The same applies to workers effect and workers trait. These two effects (=traits) will give rise to honey yield trait. We will assume that colony honey yield is approximately normally distributed with the mean of 20 kg and variance of 4 \(kg^2\), which implies that most colonies will have honey yield between 14 kg and 26 kg (see hist(rnorm(n = 1000, mean = 20, sd = sqrt(4)))). Traits like honey yield have a complex polygenic genetic architecture, so we will assume that this trait is influenced by 100 QTL per chromosome (with 16 chromosomes, this gives us 1600 QTL in total).

We will first initiate global simulation parameters and set the mean of queen effects to 10 kg with genetic variance of 1 \(kg^2\), while we will set the mean of workers effects to 10 kg with genetic variance of 1 \(kg^2\). The mean and the variance for the worker effect are proportionally scaled by the expected number of workers in a colony. The mean and variance for the queen effect is assumed larger than for the workers effect, because there is one queen and many workers in colony and we assume that workers effects “accumulate”. Deciding how to split the colony mean between queen and workers effects will depend on the individual to colony mapping function, which we will describe in the Colony value sub-section.

# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

nQtlPerChr <- 100

# Genetic parameters for queen and workers effects - each represented by a trait
mean <- c(10, 10 / SP$nWorkers)
varA <- c(1, 1 / SP$nWorkers)

We next set genetic correlation between the queen and workers effects to -0.5 to reflect the commonly observed antagonistic relationship between these effects. With all the quantitative genetic parameters defined, we now add two additive traits to global simulation parameters and name them queenTrait and workerTrait. These parameters drive the simulation of QTL effects. Read about all the other trait simulation options in AlphaSimR via: vignette(topic = "traits", package="AlphaSimR").

corA <- matrix(data = c( 1.0, -0.5, 
                        -0.5,  1.0), nrow = 2, byrow = TRUE)
SP$addTraitA(nQtlPerChr = nQtlPerChr, mean = mean, var = varA, corA = corA,
             name = c("queenTrait", "workersTrait"))

Finally, we set the environmental variance of the queen and workers effects to 3 \(kg^2\) and we again scale the worker variance by the expected number of workers. Contrary to the negative genetic correlation, we here assume that environmental correlation between the queen and workers effects is slightly positive, 0.3. This is just an example! These parameters should be based on literature or simulation scenarios of interest.

varE <- c(3, 3 / SP$nWorkers)
corE <- matrix(data = c(1.0, 0.3, 
                        0.3, 1.0), nrow = 2, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)

Individual values in the base population

Now we create a base population of virgin queens. Since we defined two traits, all honeybees in the simulation will have genetic and phenotypic values for both traits. The genetic values are stored in the gv slot of each Pop object, while phenotypic values are stored in the pheno slot.

#>      queenTrait workersTrait
#> [1,]   9.223527   0.11691699
#> [2,]   9.887075   0.02692530
#> [3,]   7.485241   0.20623985
#> [4,]  11.409446   0.04572082
#> [5,]   9.541581  -0.03860374
#> [6,]  10.825529   0.10350462
#>      queenTrait workersTrait
#> [1,]   6.772158   -0.1776047
#> [2,]  11.173154    0.1904641
#> [3,]   5.619330   -0.1267459
#> [4,]   9.086497   -0.1873403
#> [5,]   7.353185   -0.2317826
#> [6,]  10.278738    0.2180535

Note that these are virgin queens, yet we obtained queen and workers effect values for them! Is this wrong? No! Virgin queens carry DNA with genes that are differentially expressed in different castes, which would be only showed in their phenotype. Hence, virgin queens have genetic values for the queen and worker effects, but they might never actually express these effects. In this simulation virgin queens also obtained phenotypic values for both of the effects. This is technically incorrect because virgin queens don’t express genes for the worker effect at all, and they also do not express the queen effect, not until they become the queen of a colony. We can treat these phenotypic values for virgin queens as values that we could see if these virgin queens would express these traits. We will show later in the Colony value sub-section how we use these traits from different castes. If existence of these phenotypic values for certain castes is a hindrance, we can always remove them for population or colony objects by modifying the corresponding slots as required.

As with the virgin queens, drones also carry DNA with genes that are expressed in different castes. Therefore, drones will also have the queen and workers effect genetic (and phenotypic values) for honey yield even though they do not contribute to this trait in a colony.

#>      queenTrait workersTrait
#> [1,]   7.809303   0.12929587
#> [2,]   7.855408   0.11775230
#> [3,]   9.939779   0.14653540
#> [4,]   8.844980   0.11467852
#> [5,]  10.521163  -0.10273822
#> [6,]  10.364750  -0.08528203

Individual values in a colony

We continue by creating a colony from one base population virgin queen, crossing it, and adding some workers.

colony <- createColony(x = basePop[6])
colony <- cross(x = colony, drones = drones, checkCross = "warning")
colony <- addWorkers(x = colony, nInd = 50)
colony
#> An object of class "Colony" 
#> Id: 1 
#> Location: 0 0 
#> Queen: 6 
#> Number of fathers: 15 
#> Number of workers: 50 
#> Number of drones: 0 
#> Number of virgin queens: 0 
#> Has split: FALSE 
#> Has swarmed: FALSE 
#> Has superseded: FALSE 
#> Has collapsed: FALSE 
#> Is productive: FALSE

We can access the genetic and phenotypic values of colony members with functions getGv() and getPheno(), both of which have the caste argument (see more via help(getGv)).

getGv(colony, caste = "queen")
#>   queenTrait workersTrait
#> 6   10.82553    0.1035046
getGv(colony, caste = "workers") |> head(n = 4)
#>    queenTrait workersTrait
#> 36   9.047101   0.19024810
#> 37   8.772361   0.09361389
#> 38   9.435015   0.15220995
#> 39  11.040088  -0.03436708

getPheno(colony, caste = "queen")
#>   queenTrait workersTrait
#> 6   10.27874    0.2180535
getPheno(colony, caste = "workers") |> head(n = 4)
#>    queenTrait workersTrait
#> 36   8.586775   0.12283613
#> 37   8.932267  -0.08197344
#> 38  11.786789   0.29553087
#> 39  10.427982  -0.06841102

For convenience, there are also alias functions for accessing the genetic and phenotypic values of each caste directly.

getQueenGv(colony)
#>   queenTrait workersTrait
#> 6   10.82553    0.1035046
getWorkersGv(colony) |> head(n = 4)
#>    queenTrait workersTrait
#> 36   9.047101   0.19024810
#> 37   8.772361   0.09361389
#> 38   9.435015   0.15220995
#> 39  11.040088  -0.03436708

getQueenPheno(colony)
#>   queenTrait workersTrait
#> 6   10.27874    0.2180535
getWorkersPheno(colony) |> head(n = 4)
#>    queenTrait workersTrait
#> 36   8.586775   0.12283613
#> 37   8.932267  -0.08197344
#> 38  11.786789   0.29553087
#> 39  10.427982  -0.06841102

Some phenotypes, such as honey yield, are only expressed if colony is at full size. This is achieved by the buildUp() colony event function that adds worker and drones and hence turns on the production status of the colony (to TRUE). SIMplyBee includes a function ìsProductive() to check the production status of a colony.

# Check if colony is productive
isProductive(colony)
#> [1] FALSE

# Build-up the colony and check the production status again
colony <- buildUp(colony)
colony
#> An object of class "Colony" 
#> Id: 1 
#> Location: 0 0 
#> Queen: 6 
#> Number of fathers: 15 
#> Number of workers: 100 
#> Number of drones: 100 
#> Number of virgin queens: 0 
#> Has split: FALSE 
#> Has swarmed: FALSE 
#> Has superseded: FALSE 
#> Has collapsed: FALSE 
#> Is productive: TRUE
isProductive(colony)
#> [1] TRUE

For the ease of further demonstration, we now combine workers’ values into a single data.frame.

# Collate genetic and phenotypic values of workers
df <- data.frame(id = colony@workers@id,
                 mother = colony@workers@mother,
                 father = colony@workers@father,
                 gvQueenTrait = colony@workers@gv[, "queenTrait"],
                 gvWorkersTrait = colony@workers@gv[, "workersTrait"],
                 pvQueenTrait =  colony@workers@pheno[, "queenTrait"],
                 pvWorkersTrait = colony@workers@pheno[, "workersTrait"])
head(df)
#>   id mother father gvQueenTrait gvWorkersTrait pvQueenTrait pvWorkersTrait
#> 1 86      6     21     8.614909     0.21184920     7.143641    0.236943661
#> 2 87      6     24     9.471665     0.12572729    11.683815    0.354132710
#> 3 88      6     33     9.320402     0.12789536    10.291714   -0.026447705
#> 4 89      6     25    10.247971    -0.01215820     7.563013   -0.002639673
#> 5 90      6     23    10.254895     0.26071445     9.431241    0.211378922
#> 6 91      6     22     9.414942     0.08933521     7.958763   -0.095667911

To visualise correlation between queen and workers effects in workers, we plot these effect values against each other.

# Covariation between queen and workers effect genetic values in workers
p <- ggplot(data = df, aes(x = gvQueenTrait, y = gvWorkersTrait)) +
  xlab("Genetic value for the queen effect") +
  ylab("Genetic value for the workers effect") +
  geom_point() +
  theme_classic()
print(p)

In SIMplyBee, we know genetic values of all individuals, including drones that the queen mated with (=fathers in a colony)!

# Variation in patriline genetic values
getFathersGv(colony)
#>    queenTrait workersTrait
#> 21   7.809303  0.129295867
#> 22   7.855408  0.117752302
#> 23   9.939779  0.146535404
#> 24   8.844980  0.114678515
#> 25  10.521163 -0.102738217
#> 26  10.364750 -0.085282028
#> 27   7.646466  0.152260317
#> 28   7.403455  0.239805424
#> 29   6.770573  0.225434636
#> 30  12.680537 -0.043309584
#> 31  11.387844 -0.023457182
#> 32  11.998622  0.003557936
#> 33   8.020824  0.094440882
#> 34   9.376001  0.028725672
#> 35   8.782735  0.073351749

Knowing the father of each worker, we inspect variation in the distribution of genetic values of worker by the patriline (workers from a single father drone) for the workers effect.

Colony value

However, in honeybees we usually don’t observe values on individuals, but on a colony. SIMplyBee provides functions for mapping individual values to a colony value. The general function for this is calcColonyValue(), which can combine any value and trait from any caste. There are also aliases calcColonyGv() and calcColonyPheno(). These functions require users to specify the so-called mapping function (via the FUN argument). The mapping function specifies queen and workers traits (potentially also drone traits) and what function we want to apply to each of them before mapping them to the colony value(s). We can also specify whether the colony value(s) depend on the production status. For example, if a colony is not productive, its honey yield would be 0 or unobserved. SIMplyBee provides a general mapping function mapCasteToColonyValue() and aliases mapCasteToColonyGv() and mapCasteToColonyPheno(). These functions have arguments to cater for various situations. By default, they first calculate caste values: leave the queen’s value as it is, sum workers’ values, potentially sum drones’ values, and lastly sum all these caste values together into a colony value. Users can provide their own mapping function(s) too!

We now calculate honey yield for our colony - a single value for the colony.

# Colony phenotype value
calcColonyPheno(colony, queenTrait = "queenTrait", workersTrait = "workersTrait")
#>          [,1]
#> [1,] 16.89724
help(calcColonyPheno)
help(mapCasteToColonyPheno)

These colony values are not stored in a colony, because they change as colony changes due to various events. For example, reducing the number of workers will reduce the colony honey yield.

# Colony phenotype value from a reduced colony
removeWorkers(colony, p = 0.5) |>
  calcColonyPheno(queenTrait = "queenTrait", workersTrait = "workersTrait")
#>          [,1]
#> [1,] 13.51124

Please note that we assumed that the queen contributes half to colony honey yield and workers contribute the other half. This means that removing workers will still give a non-zero honey yield! This shows that we have to design the mapping between individual, caste, and colony values with care!

# Colony phenotype value from a reduced colony
removeWorkers(colony, p = 0.99) |>
  calcColonyPheno(queenTrait = "queenTrait", workersTrait = "workersTrait")
#>          [,1]
#> [1,] 10.10042

Finally, note that SIMplyBee currently does not provide functionality for breeding values, dominance deviations, and epistatic deviations at caste and colony levels, despite the availabiliy of AlphaSimR bv(), dd(), and aa() functions. This is because we have to check or develop theory on how to calculate these values across active colonies and hence we currently advise against the use of AlphaSimR bv(), dd(), and aa() functions with SIMplyBee as the output of these functions could be easily misinterpreted.

MultiColony values

The same functions can be used on a MultiColony class object. Let’s create an apiary.

apiary <- createMultiColony(basePop[7:20])
drones <- createDrones(basePop[1:5], nInd = 100)
droneGroups <- pullDroneGroupsFromDCA(drones, n = nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")
apiary <- buildUp(apiary)

We can extract the genetic and phenotypic values from multiple colonies in the same manner as from a single colony, by using get*Gv() and get*Pheno() functions. The output of these function is a named list with values for each colony or a single matrix if we set the collapse argument to TRUE.

getQueenGv(apiary) |> head(n = 4)
#> $`2`
#>   queenTrait workersTrait
#> 7   8.988466  -0.08005247
#> 
#> $`3`
#>   queenTrait workersTrait
#> 8   10.44032    0.1407429
#> 
#> $`4`
#>   queenTrait workersTrait
#> 9   9.434017    0.1310034
#> 
#> $`5`
#>    queenTrait workersTrait
#> 10   12.19888   0.09906585
getQueenGv(apiary, collapse = TRUE) |> head(n = 4)
#>    queenTrait workersTrait
#> 7    8.988466  -0.08005247
#> 8   10.440317   0.14074290
#> 9    9.434017   0.13100337
#> 10  12.198880   0.09906585

In a similar manner, we can calculate colony value for all the colonies in our apiary, where the row names of the output represent colony IDs.

colonyGv <- calcColonyGv(apiary)
colonyPheno <- calcColonyPheno(apiary)
data.frame(colonyGv, colonyPheno)
#>     colonyGv colonyPheno
#> 2   8.802483   10.346684
#> 3  19.783259   20.698010
#> 4  20.241753   23.257388
#> 5  19.935652   19.247002
#> 6  22.106762   21.378173
#> 7  24.147570   25.338415
#> 8  16.094285   19.183372
#> 9  14.438280   11.439974
#> 10  8.366858    2.963799
#> 11 26.407677   22.973798
#> 12 13.988678   16.177346
#> 13 20.850452   17.886973
#> 14 19.808384   18.207366
#> 15 25.280465   26.467217

Selection on colony values

Since the aim of selection is to select the best individuals or colonies for the reproduction, we could select the best colony in our apiary based on either genetic or phenotypic value for grafting the new generation of virgin queens. We can use the function selectColonies() that takes a matrix of colony values (the output of calcColonyValue() function). The default behavior is to select the colonies with the highest value (argument selectTop set to TRUE), but you can also select the colonies with the lowest values (argument selectTop set to FALSE).

# Select the best colony based on gv
selectColonies(apiary, n = 1, by = colonyGv)
#> An object of class "MultiColony" 
#> Number of colonies: 1 
#> Are empty: 0 
#> Are NULL: 0 
#> Have split: 0 
#> Have swarmed: 0 
#> Have superseded: 0 
#> Have collapsed: 0 
#> Are productive: 0
# Select the best colony based on phenotype
selectColonies(apiary, n = 1, by = colonyPheno)
#> An object of class "MultiColony" 
#> Number of colonies: 1 
#> Are empty: 0 
#> Are NULL: 0 
#> Have split: 0 
#> Have swarmed: 0 
#> Have superseded: 0 
#> Have collapsed: 0 
#> Are productive: 0

The same functionality is implemented in pullColonies() and removeColonies().

Honey yield and Calmness

In this section we expand simulation to two uncorrelated colony traits with queen and workers effects, honey yield and calmness. We follow the same recipe as in the previous section where we simulated only one colony trait.

We first reinitialize the global simulation parameters because we will define new traits. For honey yield we will use the same parameters as before, while for calmness trait we will assume that the trait is scored continuously in such a way that negative values are undesirable and positive values are desirable with zero being population mean. We will further assume the same variances for calmness as for honey yield, and a genetic (and environmental) correlation between the queen and workers effects of -0.4 (and 0.2) for calmness. We assume no genetic or environmental correlation between honey yield and calmness. Beware, this is just an example to show you how to simulate multiple colony traits - we have made up these parameters - please use literature estimates in your simulations!

# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

nQtlPerChr <- 100

# Quantitative genetic parameters - for two traits, each with the queen and workers effects
meanP <- c(10, 10 / SP$nWorkers, 0, 0)
varA <- c(1, 1 / SP$nWorkers, 1, 1 / SP$nWorkers)
corA <- matrix(data = c( 1.0, -0.5,  0.0,  0.0, 
                        -0.5,  1.0,  0.0,  0.0,
                         0.0,  0.0,  1.0, -0.4, 
                         0.0,  0.0, -0.4,  1.0), nrow = 4, byrow = TRUE)
SP$addTraitA(nQtlPerChr = 100, mean = meanP, var = varA, corA = corA,
             name = c("yieldQueenTrait", "yieldWorkersTrait",
                      "calmQueenTrait", "calmWorkersTrait"))

varE <- c(3, 3 / SP$nWorkers, 3, 3 / SP$nWorkers)
corE <- matrix(data = c(1.0, 0.3, 0.0, 0.0,
                        0.3, 1.0, 0.0, 0.0,
                        0.0, 0.0, 1.0, 0.2,
                        0.0, 0.0, 0.2, 1.0), nrow = 4, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)

We continue by creating a base population of virgin queens and from them an apiary with 10 full-sized colonies.

basePop <- createVirginQueens(founderGenomes)
drones <- createDrones(x = basePop[1:5], nInd = 100)
apiary <- createMultiColony(basePop[6:20])
droneGroups <- pullDroneGroupsFromDCA(drones, nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")
apiary <- buildUp(apiary)
apiary
#> An object of class "MultiColony" 
#> Number of colonies: 15 
#> Are empty: 0 
#> Are NULL: 0 
#> Have split: 0 
#> Have swarmed: 0 
#> Have superseded: 0 
#> Have collapsed: 0 
#> Are productive: 15

We can again inspect the genetic (and phenotypic) values of all individuals in each colony and whole apiary with get*Gv() and get*Pheno() functions. Now, the output contains four traits representing the queen and workers effect for honey yield and calmness. These functions also take an nInd argument to sample a number of individuals along with their values.

getQueenGv(apiary) |> head(n = 4)
#> $`1`
#>   yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 6        9.038413         0.1342118    0.009992128       -0.1053888
#> 
#> $`2`
#>   yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 7        10.81587       0.004706908      0.2483213        0.1355696
#> 
#> $`3`
#>   yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 8        10.19955        0.02775753     -0.8040684      -0.02174428
#> 
#> $`4`
#>   yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 9         7.39489         0.2588803     -0.1823955       -0.1214185
getWorkersPheno(apiary, nInd = 3) |> head(n = 4)
#> $`1`
#>     yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 521        9.273155       0.005411845     0.35273507       0.01918771
#> 522        9.614122       0.134026230     0.04693592       0.13083410
#> 523       11.683053      -0.049007193    -1.22413135      -0.42308038
#> 
#> $`2`
#>     yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 721       10.350785       -0.05076961       2.645384       0.49622351
#> 722        9.882067        0.13289452       1.239964       0.08641623
#> 723        7.982978       -0.35365089       2.553985       0.25026644
#> 
#> $`3`
#>     yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 921        12.18280        0.24670995      0.4901567      -0.06073114
#> 922        11.75166       -0.05321984     -1.1476476      -0.05287263
#> 923        12.01447        0.31518766     -1.7783011      -0.01929232
#> 
#> $`4`
#>      yieldQueenTrait yieldWorkersTrait calmQueenTrait calmWorkersTrait
#> 1121       11.492588         0.2318531     -0.9164149       0.05250092
#> 1122        7.426122         0.3122308     -1.4047680      -0.11462932
#> 1123        5.944146         0.2463056      1.4131271      -0.55854670

Now, we calculate colony genetic and phenotypic values for all colonies in the apiary. Since we are simulating two traits, honey yield and calmness, we have two ways to calculate corresponding colony values. The first way is to use the default mapCasteToColony*() function in calcColony*() and only define additional arguments as shown here:

colonyValues <- calcColonyPheno(apiary,
                                queenTrait = c("yieldQueenTrait", "calmQueenTrait"),
                                workersTrait = c("yieldWorkersTrait", "calmWorkersTrait"),
                                traitName = c("yield", "calmness"),
                                checkProduction = c(TRUE, FALSE)) |> as.data.frame()
colonyValues
#>       yield   calmness
#> 1  19.11033 -5.4034262
#> 2  16.13080  5.4381313
#> 3  18.32804  0.8157093
#> 4  27.23456 -4.8212366
#> 5  21.97054 -9.5186342
#> 6  30.66659 -6.3438602
#> 7  18.54297  1.7312849
#> 8  21.23205 -0.3518961
#> 9  24.41140 -5.1047658
#> 10 23.38081  2.5662358
#> 11 29.83958  0.5917187
#> 12 11.60884  3.6631743
#> 13 24.61257  5.6376562
#> 14 17.92309 -6.1450422
#> 15 12.74938 -0.1792312

The second way is to create our own mapping function. An equivalent outcome to the above is shown below just to demonstrate use of your own function, but we are simply just reusing mapCasteToColonyPheno() twice;)

myMapCasteToColonyPheno <- function(colony) {
  yield <- mapCasteToColonyPheno(colony,
                                 queenTrait = "yieldQueenTrait",
                                 workersTrait = "yieldWorkersTrait",
                                 traitName = "yield",
                                 checkProduction = TRUE)
  calmness <- mapCasteToColonyPheno(colony,
                                    queenTrait = "calmQueenTrait",
                                    workersTrait = "calmWorkersTrait",
                                    traitName = "calmness",
                                    checkProduction = FALSE)
  return(cbind(yield, calmness))
}
colonyValues <- calcColonyPheno(apiary, FUN = myMapCasteToColonyPheno) |> as.data.frame()
colonyValues
#>       yield   calmness
#> 1  19.11033 -5.4034262
#> 2  16.13080  5.4381313
#> 3  18.32804  0.8157093
#> 4  27.23456 -4.8212366
#> 5  21.97054 -9.5186342
#> 6  30.66659 -6.3438602
#> 7  18.54297  1.7312849
#> 8  21.23205 -0.3518961
#> 9  24.41140 -5.1047658
#> 10 23.38081  2.5662358
#> 11 29.83958  0.5917187
#> 12 11.60884  3.6631743
#> 13 24.61257  5.6376562
#> 14 17.92309 -6.1450422
#> 15 12.74938 -0.1792312

Again, we can now select the best colony based on the best phenotypic value for either yield, calmness, or an index of both. Let’s say that both traits are equally important so we select on a weighted sum of both of them - we will use the AlphaSimR selIndex() function that enables this calculation along with scaling. We will represent the index such that it has a mean of 100 and standard deviation of 10 units.

colonyValues$Index <- selIndex(Y = colonyValues, b = c(0.5, 0.5), scale = TRUE) * 10 + 100
bestColony <- selectColonies(apiary, n = 1, by = colonyValues$Index)
getId(bestColony)
#> [1] 13

We see that we selected colony with ID “4”, but we would be selecting a different colony based on different selection criteria (yield, calmness, or index).

Strength and honey yield

In this section we change simulation to two traits where the phenotype realisation of the first trait affects the phenotype realisation of the second trait. Specifically, we will assume that queen’s fecundity, and hence the number of workers, is under the genetic affect of the queen and her environment. Furthermore, we will assume as before that colony honey yield is due to the queen effect and workers effect. Since the value of the workers effect depends on then number of workers, we obtain correlation between fecundity and honey yield, even if these traits would be uncorrelated on the queen level. We emphasise that this is just an example and the biology of these traits might be different.

We follow the same logic as before and simulate three traits that will contribute to two colony traits, queen’s fecundity, that is colony strength, and honey yield. We assume that fecundity is only due to the queen (and not the workers), hence we simulate only the queen effect for this trait. For honey yield we again assume that both the queen and workers contribute to the colony value. For speed of simulation we only simulate 100 workers per colony on average and split honey yield mean between the queen and workers. We measure fecundity with the number of workers, which is a count variable and for such variables Poisson distribution is a good model. This distribution has just one parameter (lambda) that represents both the mean and variance of the variable. To this end we set phenotypic variance to 100 and split it into 25 for genetic and 65 for environmental variance. As before we warn that these are just exemplary values to demonstrate the code functionality and do not necessarily reflect published values!

# Global simulation parameters
SP <- SimParamBee$new(founderGenomes)

# Quantitative genetic parameters
# - the first trait has only the queen effect
# - the second trait has both the queen and workers effects
nWorkers <- 100
mean <- c(nWorkers, 10, 10 / nWorkers)
varA <- c(25, 1, 1 / nWorkers)
corA <- matrix(data = c(1.0,  0.0,  0.0,
                        0.0,  1.0, -0.5, 
                        0.0, -0.5,  1.0), nrow = 3, byrow = TRUE)
SP$addTraitA(nQtlPerChr = 100, mean = mean, var = varA, corA = corA,
             name = c("fecundityQueenTrait", "yieldQueenTrait", "yieldWorkersTrait"))

varE <- c(75, 3, 3 / nWorkers)
corE <- matrix(data = c(1.0, 0.0, 0.0,
                        0.0, 1.0, 0.3,
                        0.0, 0.3, 1.0), nrow = 3, byrow = TRUE)
SP$setVarE(varE = varE, corE = corE)

We continue by creating an apiary with 10 colonies.

basePop <- createVirginQueens(founderGenomes)
drones <- createDrones(x = basePop[1:5], nInd = 100)
apiary <- createMultiColony(basePop[6:20])
droneGroups <- pullDroneGroupsFromDCA(drones, nColonies(apiary), nDrones = 15)
apiary <- cross(x = apiary, drones = droneGroups, checkCross = "warning")

Let’s explore queen’s genetic and phenotypic values for fecundity and honey yield. The below printouts show quite some variation in fecundity between queens at the genetic, but particularly phenotypic level. This is a small example, so we should not put too much into correlations between these three variables. However, if you restart this simulation many times, you will notice zero correlation on average between fecundityQueenTrait and the other two traits and negative correlation on average between yieldQueenTrait and yieldWorkersTrait. Just like we defined in the global simulation parameters.

#>    fecundityQueenTrait yieldQueenTrait yieldWorkersTrait
#> 6            107.75417        8.654938       0.344653757
#> 7            101.78639       12.412649      -0.034063855
#> 8             99.62434       10.281212       0.222244847
#> 9             95.29052       10.000812      -0.055682404
#> 10            95.23090       10.807911       0.079380043
#> 11            91.14948       10.385416      -0.003545913
#> 12           101.79576        9.123699       0.115922378
#> 13            91.70656       10.722697       0.080830618
#> 14           111.35943       10.508601      -0.011382206
#> 15            99.88194        9.419631       0.143265980
#> 16           100.46803        8.455684       0.037671856
#> 17            96.71180       10.482327       0.106537139
#> 18            96.35590       10.619149       0.034493450
#> 19           100.66377       11.544070       0.062911576
#> 20            96.23236       10.130749       0.069959028
#>                     fecundityQueenTrait yieldQueenTrait yieldWorkersTrait
#> fecundityQueenTrait           1.0000000       0.4025526        -0.0385327
#> yieldQueenTrait               0.4025526       1.0000000         0.0990734
#> yieldWorkersTrait            -0.0385327       0.0990734         1.0000000

We next build-up colonies in the apiary. But instead of building them all up to the same fixed number of workers, we build them up according to queen’s fecundity. For that we use the sampling function nWorkersColonyPhenotype(), that samples the number of workers based on phenotypes of colony members, in our case fecundityQueenTrait in queens. Correspondingly, each colony will have a different number of workers. Read more about this function in it’s help page.

apiary <- buildUp(apiary, nWorkers = nWorkersColonyPhenotype,
                  queenTrait = "fecundityQueenTrait")
cbind(nWorkers = nWorkers(apiary), queenPheno)
#>    nWorkers fecundityQueenTrait yieldQueenTrait yieldWorkersTrait
#> 1        98            97.75481        8.375729        0.25512271
#> 2       100           108.20195       15.317456       -0.18684285
#> 3        92            91.82220        9.241177        0.17633134
#> 4       103           102.86006        9.167062       -0.21767959
#> 5        89            89.37903        9.616343       -0.13231939
#> 6        92            92.11085       12.204219        0.11633735
#> 7       102           101.64082        9.999383       -0.07570675
#> 8        98            97.62292       10.732032        0.38565160
#> 9       106           106.09796       10.571290       -0.17443362
#> 10       99            99.46442       13.094534        0.54966329
#> 11       86            85.60070        9.153848       -0.04867843
#> 12       89            89.25824        9.850446       -0.21275894
#> 13       97            97.43110        9.582485       -0.22666368
#> 14      101           101.04621        9.004872        0.12743674
#> 15      100            99.74655        9.960594        0.25673993
help(nWorkersColonyPhenotype)

To compute the colony value for honey yield, we again employ the calcColonyPheno() function. Correlating the queen and colony values we will now see a positive correlation because our individual to colony mapping function sums workers effect across all workers and the more workers there are the larger the sum.

#>                        nWorkers fecundityQueenTrait yieldQueenTrait
#> nWorkers             1.00000000           0.9452944      0.16471597
#> fecundityQueenTrait  0.94529435           1.0000000      0.40255260
#> yieldQueenTrait      0.16471597           0.4025526      1.00000000
#> yieldWorkersTrait    0.05707226          -0.0385327      0.09907340
#> yield               -0.19774614          -0.2445052     -0.08927866
#>                     yieldWorkersTrait       yield
#> nWorkers                   0.05707226 -0.19774614
#> fecundityQueenTrait       -0.03853270 -0.24450518
#> yieldQueenTrait            0.09907340 -0.08927866
#> yieldWorkersTrait          1.00000000  0.60039034
#> yield                      0.60039034  1.00000000

Quantitative genetics

2024-09-19

Introduction