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:

1. the number of traits,
2. the number of quantitative trait loci (QTL) that affect the traits,
3. the distribution of QTL effects,
4. trait means, and
5. 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², 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², while we will set the mean of workers effects to 10 kg with genetic variance of 1 kg². 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² 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,]  11.491123  0.010006062
#> [2,]   9.158856  0.039861300
#> [3,]  11.178754  0.001518615
#> [4,]  10.253542  0.070062259
#> [5,]   9.102163  0.031502155
#> [6,]   9.610434  0.102724424
#>      queenTrait workersTrait
#> [1,]  11.444570   0.15109914
#> [2,]   9.689239   0.13899686
#> [3,]  11.929311  -0.20162124
#> [4,]  11.116019   0.11314917
#> [5,]   6.417610   0.11798205
#> [6,]   8.566441   0.07864181

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,]  13.463570   0.10184425
#> [2,]  11.081738   0.08776253
#> [3,]  13.019753  -0.13721622
#> [4,]   7.644546   0.02320137
#> [5,]   9.838458  -0.03455270
#> [6,]   8.285439   0.08863234

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   9.610434    0.1027244
getGv(colony, caste = "workers") |> head(n = 4)
#>    queenTrait workersTrait
#> 36   9.850256  -0.01993718
#> 37  11.411114   0.05160111
#> 38   8.769486   0.06829881
#> 39   7.962571   0.09183652

getPheno(colony, caste = "queen")
#>   queenTrait workersTrait
#> 6   8.566441   0.07864181
getPheno(colony, caste = "workers") |> head(n = 4)
#>    queenTrait workersTrait
#> 36   8.210533  -0.20002026
#> 37  10.690111  -0.28236479
#> 38   7.785036  -0.14989946
#> 39   6.967794   0.05981064

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

getQueenGv(colony)
#>   queenTrait workersTrait
#> 6   9.610434    0.1027244
getWorkersGv(colony) |> head(n = 4)
#>    queenTrait workersTrait
#> 36   9.850256  -0.01993718
#> 37  11.411114   0.05160111
#> 38   8.769486   0.06829881
#> 39   7.962571   0.09183652

getQueenPheno(colony)
#>   queenTrait workersTrait
#> 6   8.566441   0.07864181
getWorkersPheno(colony) |> head(n = 4)
#>    queenTrait workersTrait
#> 36   8.210533  -0.20002026
#> 37  10.690111  -0.28236479
#> 38   7.785036  -0.14989946
#> 39   6.967794   0.05981064

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     35     9.630311     0.07334735     8.425254     0.13203156
#> 2 87      6     34     9.725538    -0.03805660    11.150475    -0.16073869
#> 3 88      6     28     9.338029     0.05223589     7.671042    -0.03789632
#> 4 89      6     31    10.022829     0.01472072     9.652027     0.05032514
#> 5 90      6     29    10.842829     0.05126267    13.379729     0.40398203
#> 6 91      6     26     9.136576     0.07576119     7.637438    -0.09484745

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  13.463570  0.101844254
#> 22  11.081738  0.087762529
#> 23  13.019753 -0.137216220
#> 24   7.644546  0.023201372
#> 25   9.838458 -0.034552704
#> 26   8.285439  0.088632338
#> 27  10.355917  0.123242748
#> 28   9.619326 -0.063217229
#> 29  12.653896 -0.021486725
#> 30  10.284908 -0.099095455
#> 31   9.503393  0.057405645
#> 32  11.310419  0.124040629
#> 33   9.082032  0.038273491
#> 34   9.803302 -0.001359011
#> 35   9.120424  0.080232423

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,] 17.63344
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,] 12.44763

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,] 8.657439

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   9.350515     0.161232
#> 
#> $`3`
#>   queenTrait workersTrait
#> 8   10.05018     0.131991
#> 
#> $`4`
#>   queenTrait workersTrait
#> 9   9.475807    0.1251078
#> 
#> $`5`
#>    queenTrait workersTrait
#> 10   10.92759   0.09629173
getQueenGv(apiary, collapse = TRUE) |> head(n = 4)
#>    queenTrait workersTrait
#> 7    9.350515   0.16123201
#> 8   10.050178   0.13199104
#> 9    9.475807   0.12510775
#> 10  10.927588   0.09629173

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  20.87496   21.960813
#> 3  17.69244   19.653500
#> 4  20.08088   20.996948
#> 5  17.71601   12.113224
#> 6  15.50407   11.561151
#> 7  22.94692   23.166751
#> 8  18.65563   18.460645
#> 9  10.48828   15.888038
#> 10 22.47712   27.329221
#> 11 13.03558   12.847836
#> 12  6.64320    4.680641
#> 13 24.30971   22.760854
#> 14 20.36717   22.464717
#> 15 21.07002   20.937064

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().