diff --git a/.gitignore b/.gitignore
index c8dd2d6a..7a82f371 100644
--- a/.gitignore
+++ b/.gitignore
@@ -40,3 +40,5 @@ vignettes/*.pdf
 
 # Other
 .DS_Store
+src/*.o
+src/*.so
diff --git a/DESCRIPTION b/DESCRIPTION
index 5ac2bc7d..94b58e24 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: SIMplyBee
 Type: Package
 Title: 'AlphaSimR' Extension for Simulating Honeybee Populations and Breeding Programmes
-Version: 0.4.1
+Version: 0.4.2
 Authors@R: c(
   person("Jana", "Obšteter", email = "obsteter.jana@gmail.com",
   role = c("aut", "cre"), comment = c(ORCID = "0000-0003-1511-3916")),
@@ -27,7 +27,7 @@ License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
 Imports: methods, R6, stats, utils, extraDistr (>= 1.9.1), RANN, Rcpp (>= 0.12.7)
-Depends: R (>= 3.3.0), AlphaSimR (>= 1.5.3)
+Depends: R (>= 3.3.0), AlphaSimR (>= 1.6.1)
 LinkingTo: Rcpp, RcppArmadillo (>= 0.7.500.0.0), BH
 RoxygenNote: 7.3.2
 Suggests: 
diff --git a/NAMESPACE b/NAMESPACE
index 231169a0..07ec9e03 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -69,7 +69,6 @@ export(getGv)
 export(getIbdHaplo)
 export(getId)
 export(getLocation)
-export(getMisc)
 export(getPheno)
 export(getPooledGeno)
 export(getQtlGeno)
@@ -192,7 +191,6 @@ export(replaceWorkers)
 export(resetEvents)
 export(selectColonies)
 export(setLocation)
-export(setMisc)
 export(setQueensYearOfBirth)
 export(simulateHoneyBeeGenomes)
 export(split)
diff --git a/NEWS.md b/NEWS.md
index d72a41ac..33dbee69 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,28 @@ editor_options:
     wrap: 72
 ---
 
+# SIMplyBee version 0.4.2
+
+-   20??-??-??
+
+## Major changes
+-   TODO
+
+## New features
+-   TODO
+
+## Bug fixes
+-   editCsdLocus() works now also with just 2 alleles #591
+
+## Backgrond/package/etc. work
+-   calcBeeGRMIbs() can now return centred genotype matrix, allele frequencies,
+and scale factor used to calculate the GRM #594
+-   Improved default for getPooledGeno() (to type="mean"") and added an example
+on how to collect pooled workers' genotype accross colonies (but pooling is
+done within a colony!) #592
+-   We now removed setMisc() and getMisc() because we now use the new AlphaSimR
+structure of the misc slot that is much easier to use. #584
+
 # SIMplyBee version 0.4.1
 
 -   2024-09-19
@@ -26,7 +48,7 @@ which caused an error. We now read in the locations from a csv file.
     now c(0, 0) PR#500
 
 
-## New features ##
+## New features
 -   In setLocation(MultiColony) we can set one location (numeric) or
     multiple (list or data.frame) PR#500
 -   getLocation(MultiColony) got the collapse argument
@@ -57,8 +79,6 @@ which caused an error. We now read in the locations from a csv file.
 -   Bug fix - get\*Haplo() functions were returning diploid drones when
     input was a Pop-class
 
--   
-
 # SIMplyBee version 0.3.0
 
 -   2022-12-05 First public/CRAN version of the package
diff --git a/R/Functions_L0_auxilary.R b/R/Functions_L0_auxilary.R
index 28e56689..3e774144 100644
--- a/R/Functions_L0_auxilary.R
+++ b/R/Functions_L0_auxilary.R
@@ -1079,11 +1079,7 @@ getQueenAge <- function(x, currentYear, simParamBee = NULL) {
     }
   } else if (isColony(x)) {
     if (isQueenPresent(x, simParamBee = simParamBee)) {
-      if(packageVersion("AlphaSimR") > package_version("1.5.3")){
-        ret <- currentYear - x@queen@misc$yearOfBirth[[1]]
-      }else{
-        ret <- currentYear - x@queen@misc[[1]]$yearOfBirth
-      }
+      ret <- currentYear - x@queen@misc$yearOfBirth[[1]]
     } else {
       ret <- NA
     }
@@ -1280,8 +1276,6 @@ getCasteId <- function(x, caste = "all", collapse = FALSE, simParamBee = NULL) {
 #'   vector with sex information
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getCaste}}
-#'
 #' @return when \code{x} is \code{\link[AlphaSimR]{Pop-class}} for \code{caste != "all"}
 #'  or list for \code{caste == "all"} with sex nodes named by caste;
 #'    when \code{x} is \code{\link[SIMplyBee]{Colony-class}} return is a named list of
@@ -2873,7 +2867,7 @@ nCsdAlleles <- function(x, collapse = FALSE, simParamBee = NULL) {
 #'   with haplotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getIbdHaplo}} and \code{\link[AlphaSimR]{pullIbdHaplo}}
+#' @seealso \code{\link[AlphaSimR]{pullIbdHaplo}}
 #'
 #' @return matrix with haplotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}}
 #'   and list of matrices with haplotypes when \code{x} is
@@ -3128,7 +3122,7 @@ getDronesIbdHaplo <- function(x, nInd = NULL, chr = NULL, snpChip = NULL,
 #'   with haplotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getQtlHaplo}} and \code{\link[AlphaSimR]{pullQtlHaplo}} as well as
+#' @seealso \code{\link[SIMplyBee]{getQtlGeno}} and \code{\link[AlphaSimR]{pullQtlHaplo}} as well as
 #'   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 #'
 #' @return matrix with haplotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}}
@@ -3378,7 +3372,9 @@ getDronesQtlHaplo <- function(x, nInd = NULL,
 #'   with genotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getQtlGeno}} and \code{\link[AlphaSimR]{pullQtlGeno}} as well as
+#' @seealso \code{\link[SIMplyBee]{getQtlHaplo}},
+#'   \code{\link[AlphaSimR]{pullQtlGeno}}, and
+#'   \code{\link[SIMplyBee]{getPooledGeno}}, as well as
 #'   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 #'
 #' @return matrix with genotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}} and
@@ -3628,7 +3624,7 @@ getDronesQtlGeno <- function(x, nInd = NULL,
 #'   with haplotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getSegSiteHaplo}} and \code{\link[AlphaSimR]{pullSegSiteHaplo}}
+#' @seealso \code{\link[SIMplyBee]{getSegSiteGeno}} and \code{\link[AlphaSimR]{pullSegSiteHaplo}}
 #'
 #' @return matrix with haplotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}}
 #'   and list of matrices with haplotypes when \code{x} is
@@ -3870,7 +3866,9 @@ getDronesSegSiteHaplo <- function(x, nInd = NULL,
 #'   with genotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getSegSiteGeno}} and \code{\link[AlphaSimR]{pullSegSiteGeno}}
+#' @seealso \code{\link[SIMplyBee]{getSegSiteHaplo}},
+#'   \code{\link[AlphaSimR]{pullSegSiteGeno}}, and
+#'   \code{\link[SIMplyBee]{getPooledGeno}}
 #'
 #' @return matrix with genotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}} and
 #'   list of matrices with genotypes when \code{x} is
@@ -4110,7 +4108,7 @@ getDronesSegSiteGeno <- function(x, nInd = NULL,
 #'   with haplotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getSnpHaplo}} and \code{\link[AlphaSimR]{pullSnpHaplo}}
+#' @seealso \code{\link[SIMplyBee]{getSnpGeno}} and \code{\link[AlphaSimR]{pullSnpHaplo}}
 #'
 #' @return matrix with haplotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}}
 #'   and list of matrices with haplotypes when \code{x} is
@@ -4351,7 +4349,9 @@ getDronesSnpHaplo <- function(x, nInd = NULL,
 #'   with genotypes of all the individuals
 #' @param simParamBee \code{\link[SIMplyBee]{SimParamBee}}, global simulation parameters
 #'
-#' @seealso \code{\link[SIMplyBee]{getSnpGeno}} and \code{\link[AlphaSimR]{pullSnpGeno}}
+#' @seealso \code{\link[SIMplyBee]{getSnpHaplo}},
+#'  \code{\link[AlphaSimR]{pullSnpGeno}}, and
+#'  \code{\link[SIMplyBee]{getPooledGeno}}
 #'
 #' @return matrix with genotypes when \code{x} is \code{\link[SIMplyBee]{Colony-class}} and
 #'   list of matrices with genotypes when \code{x} is
@@ -4577,14 +4577,14 @@ getDronesSnpGeno <- function(x, nInd = NULL,
 #'   genotypes to mimic genotyping of a pool of colony members.
 #'
 #' @param x matrix, true genotypes with individuals in rows and sites in columns
-#' @param type character, "mean" for average genotype or "count" for the counts
-#'   of reference and alternative alleles
+#' @param type character, \code{"mean"} for average genotype (default) or
+#'   \code{"count"} for the counts of reference and alternative alleles
 #' @param sex character, vector of "F" and "M" to denote the sex of individuals
 #'   in \code{x}
 #'
 #' @return a numeric vector with average allele dosage when \code{type = "mean"}
-#'   and a two-row matrix with the counts of reference (1st row) and
-#'   alternative (2nd row) alleles
+#'   or a two-row matrix with the counts of reference (1st row) and
+#'   alternative (2nd row) alleles when \code{type = "count"}
 #'
 #' @examples
 #' founderGenomes <- quickHaplo(nInd = 3, nChr = 1, segSites = 50)
@@ -4636,11 +4636,19 @@ getDronesSnpGeno <- function(x, nInd = NULL,
 #'
 #' # As an exercise you could repeat the above with different numbers of workers!
 #'
+#' # How to get pooled genotypes of workers across multiple colonies?
+#' tmp <- getSegSiteGeno(x = apiary, caste = "workers")
+#' (tmp2 = lapply(X = tmp, FUN = getPooledGeno)) # as a list of one row matrices
+#' t(sapply(X = tmp, FUN = getPooledGeno)) # as one matrix - option A
+#' do.call(what = rbind, args = tmp2) # as one matrix - option B
 #' @export
-getPooledGeno <- function(x, type = NULL, sex = NULL) {
+getPooledGeno <- function(x, type = "mean", sex = NULL) {
   if (!is.matrix(x)) {
     stop("Argument x must be a matrix class object!")
   }
+  if (is.null(type) | !(type %in% c("mean", "count"))) {
+    stop("Argument type must be specified as either mean or count!")
+  }
   n <- nrow(x)
   if (is.null(sex)) {
     warning("Argument sex is NULL. Assuming that all individuals are female/diploid!")
@@ -4658,8 +4666,6 @@ getPooledGeno <- function(x, type = NULL, sex = NULL) {
   } else if (type == "count") {
     ret <- rbind(nPloids - ret, ret)
     rownames(ret) <- c("0", "1")
-  } else {
-    stop("Argument type must be mean or count!")
   }
   return(ret)
 }
@@ -4680,9 +4686,16 @@ getPooledGeno <- function(x, type = NULL, sex = NULL) {
 #' @param sex character vector denoting sex for individuals with genotypes in
 #'   \code{x} - \code{"F"} for female and \code{"M"} for male
 #' @param alleleFreq numeric, vector of allele frequencies for the sites in
-#'   \code{x}; if \code{NULL}, then \code{\link[SIMplyBee]{calcBeeAlleleFreq}} is used
+#'   \code{x}; if \code{NULL}, then \code{\link[SIMplyBee]{calcBeeAlleleFreq}}
+#'   is used
+#' @param returnComponents logical, return GRM as well as the components used
+#'   to compute it (useful for GWAS by GBLUP)
 #'
-#' @return matrix of genomic relatedness coefficients
+#' @return if \code{returnComponents = FALSE} (default) return a matrix of
+#'   genomic relatedness coefficients; if \code{returnComponents = TRUE} return
+#'   a list with the GRM, centred genotype matrix, allele frequencies, and
+#'   scaling factor used to scale the crossproduct of centred genotype matrix
+#'   to get the GRM.
 #'
 #' @references Druet and Legarra (2020) Theoretical and empirical comparisons of
 #'   expected and realized relationships for the X-chromosome. Genetics
@@ -4708,12 +4721,12 @@ getPooledGeno <- function(x, type = NULL, sex = NULL) {
 #' GRM <- calcBeeGRMIbs(x = geno, sex = sex)
 #' # You can visualise this matrix with the function image() from the package 'Matrix'
 #'
-#' #Look at the diagonal at the relationship matrix
+#' # Look at the diagonal at the relationship matrix
 #' x <- diag(GRM)
 #' hist(x)
 #' summary(x)
 #'
-#' #Look at the off-diagonal at the relationship matrix
+#' # Look at the off-diagonal at the relationship matrix
 #' x <- GRM[lower.tri(x = GRM, diag = FALSE)]
 #' hist(x)
 #' summary(x)
@@ -4749,8 +4762,12 @@ getPooledGeno <- function(x, type = NULL, sex = NULL) {
 #' calcBeeGRMIbs(x = rbind(queenGeno, pooledGenoW), sex = c("F","F"))
 #' # You can now compare how this compare to relationships between the queen
 #' # individual workers!
+#'
+#' # Return components
+#' calcBeeGRMIbs(x = rbind(queenGeno, pooledGenoW), sex = c("F","F"),
+#'               returnComponents = TRUE)
 #' @export
-calcBeeGRMIbs <- function(x, sex, alleleFreq = NULL) {
+calcBeeGRMIbs <- function(x, sex, alleleFreq = NULL, returnComponents = FALSE) {
   if (!is.matrix(x)) {
     stop("Argument x must be a matrix class object!")
   }
@@ -4780,8 +4797,13 @@ calcBeeGRMIbs <- function(x, sex, alleleFreq = NULL) {
     #       This would overwrite x only once, at expense of doubling RAM
     x[, site] <- x[, site] - ploidy * alleleFreq[site]
   }
-  G <- tcrossprod(x) / (2 * sum(alleleFreq * (1 - alleleFreq)))
-  return(G)
+  scale <- 2 * sum(alleleFreq * (1 - alleleFreq))
+  G <- tcrossprod(x) / scale
+  if (returnComponents) {
+    return(list(G = G, x = x, alleleFreq = alleleFreq, scale = scale))
+  } else {
+    return(G)
+  }
 }
 
 #' @describeIn calcBeeGRMIbs Calculate allele frequencies from honeybee genotypes
@@ -6326,6 +6348,16 @@ calcColonyAa <- function(x, FUN = mapCasteToColonyAa, simParamBee = NULL, ...) {
 #'   in \code{\link[SIMplyBee]{SimParamBee}}. The two csd alleles must be different to
 #'   ensure heterozygosity at the csd locus.
 #' @param simParamBee global simulation parameters.
+#' @examples
+#' founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 10)
+#'
+#' SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 2)
+#' tmp <- createVirginQueens(founderGenomes)
+#' getCsdAlleles(tmp)
+#'
+#' SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 4)
+#' tmp <- createVirginQueens(founderGenomes)
+#' getCsdAlleles(tmp)
 #'
 #' @return Returns an object of \code{\link[AlphaSimR]{Pop-class}}
 editCsdLocus <- function(pop, alleles = NULL, simParamBee = NULL) {
@@ -6338,7 +6370,10 @@ editCsdLocus <- function(pop, alleles = NULL, simParamBee = NULL) {
     alleles <- expand.grid(as.data.frame(matrix(rep(0:1, length(csdSites)), nrow = 2, byrow = FALSE)))
     # Sample two different alleles (without replacement) for each individual
     nAlleles <- simParamBee$nCsdAlleles
-    alleles <- sapply(seq_len(pop@nInd), FUN = function(x) list(alleles[sample(nAlleles, size = 2, replace = F), ]))
+    alleles <- sapply(X = seq_len(pop@nInd),
+                      FUN = function(x) {
+                        list(alleles[sample(nAlleles, size = 2, replace = FALSE), , drop = FALSE])
+                      })
   }
 
   if (pop@nInd != length(alleles)) {
@@ -6548,93 +6583,6 @@ createCrossPlan <- function(x,
   return(crossPlan)
 }
 
-# Misc helpers
-# These functions replace the defunct functions of the same name in AlphaSimR
-
-#' @rdname setMisc
-#' @title Set miscellaneous information in a population
-#'
-#' @description Set miscellaneous information in a population
-#'
-#' @param x \code{\link[AlphaSimR]{Pop-class}}
-#' @param node character, name of the node to set within the \code{x@misc} slot
-#' @param value, value to be saved into \code{x@misc[[*]][[node]]}; length of
-#'   \code{value} should be equal to \code{nInd(x)}; if its length is 1, then
-#'   it is repeated using \code{rep} (see examples)
-#'
-#' @details A \code{NULL} in \code{value} is ignored
-#'
-#' @return \code{\link[AlphaSimR]{Pop-class}}
-#'
-#' @export
-setMisc <- function(x, node = NULL, value = NULL) {
-  if (isPop(x)) {
-    if (is.null(node)) {
-      stop("Argument node must be provided!")
-    }
-    if (is.null(value)) {
-      stop("Argument value must be provided!")
-    }
-    n <- nInd(x)
-    if (length(value) == 1 && n > 1) {
-      value <- rep(x = value, times = n)
-    }
-    if (length(value) != n) {
-      stop("Argument value must be of length 1 or nInd(x)!")
-    }
-
-    # Check current AlphaSimR version for new or legacy misc slot
-    if(packageVersion("AlphaSimR") > package_version("1.5.3")){
-      # New misc slot
-      x@misc[[node]] = value
-    }else{
-      # Legacy misc slot
-      names(value) = rep(x = node, times = n)
-      inode = match(names(x@misc[[1]]),node)
-      inode = inode[!is.na(inode)]
-      if(length(inode) == 0){
-        x@misc = sapply(seq_len(n),function(ind){
-          c(x@misc[[ind]],value[ind])
-        },simplify = FALSE)
-      }else{
-        x@misc = sapply(seq_len(n),function(ind){
-          c(x@misc[[ind]],value[ind])[-inode]
-        },simplify = FALSE)
-      }
-    }
-
-  }
-
-  return(x)
-}
-
-#' @rdname getMisc
-#' @title Get miscellaneous information in a population
-#'
-#' @description Get miscellaneous information in a population
-#'
-#' @param x \code{\link[AlphaSimR]{Pop-class}}
-#' @param node character, name of the node to get from the \code{x@misc} slot;
-#'   if \code{NULL} the whole \code{x@misc} slot is returned
-#'
-#' @return The \code{x@misc} slot or its nodes \code{x@misc[[*]][[node]]}
-#'
-#' @export
-getMisc <- function(x, node = NULL) {
-  if (isPop(x)) {
-    if (is.null(node)) {
-      ret <- x@misc
-    } else {
-      # Check current AlphaSimR version for new or legacy misc slot
-      ret = x@misc[[node]]
-    }
-  } else {
-    stop("Argument x must be a Pop class object!")
-  }
-  return(ret)
-}
-
-
 #' @rdname mapLoci
 #' @title Finds loci on a genetic map and return a list of positions
 #'
diff --git a/R/Functions_L1_Pop.R b/R/Functions_L1_Pop.R
index 005fe4fd..c5b823ce 100644
--- a/R/Functions_L1_Pop.R
+++ b/R/Functions_L1_Pop.R
@@ -423,10 +423,8 @@ createCastePop <- function(x, caste = NULL, nInd = NULL,
       stop("MapPop-class can only be used to create virgin queens!")
     }
     ret <- newPop(x, simParam = simParamBee)
-    if (!is.null(simParamBee$csdChr)) {
-      if (editCsd) {
-        ret <- editCsdLocus(ret, alleles = csdAlleles, simParamBee = simParamBee)
-      }
+    if (!is.null(simParamBee$csdChr) && editCsd) {
+      ret <- editCsdLocus(pop = ret, alleles = csdAlleles, simParamBee = simParamBee)
     }
     ret@sex[] <- "F"
     simParamBee$changeCaste(id = ret@id, caste = "virginQueens")
@@ -1470,15 +1468,15 @@ cross <- function(x,
         simParamBee$changeCaste(id = virginQueen@id, caste = "queen")
         simParamBee$changeCaste(id = virginQueenDrones@id, caste = "fathers")
 
-        virginQueen <- setMisc(x = virginQueen, node = "nWorkers", value = 0)
-        virginQueen <- setMisc(x = virginQueen, node = "nDrones", value = 0)
-        virginQueen <- setMisc(x = virginQueen, node = "nHomBrood", value = 0)
+        virginQueen@misc$nWorkers <- 0
+        virginQueen@misc$nDrones <- 0
+        virginQueen@misc$nHomBrood <- 0
         if (isCsdActive(simParamBee = simParamBee)) {
           val <- calcQueensPHomBrood(x = virginQueen, simParamBee = simParamBee)
         } else {
           val <- NA
         }
-        virginQueen <- setMisc(x = virginQueen, node = "pHomBrood", value = val)
+        virginQueen@misc$pHomBrood <- val
 
         if (isPop(x)) {
           ret[[virgin]] <- virginQueen
@@ -1580,20 +1578,17 @@ setQueensYearOfBirth <- function(x, year, simParamBee = NULL) {
       stop("Individuals in x must be virgin queens or queens!")
     }
     nInd <- nInd(x)
-    x <- setMisc(x = x, node = "yearOfBirth", value = year)
+    x@misc$yearOfBirth <- year
   } else if (isColony(x)) {
     if (isQueenPresent(x, simParamBee = simParamBee)) {
-      x@queen <- setMisc(x = x@queen, node = "yearOfBirth", value = year)
+      x@queen@misc$yearOfBirth <- year
     } else {
       stop("Missing queen!")
     }
   } else if (isMultiColony(x)) {
     nCol <- nColonies(x)
     for (colony in seq_len(nCol)) {
-      x[[colony]]@queen <- setMisc(
-        x = x[[colony]]@queen, node = "yearOfBirth",
-        value = year
-      )
+      x[[colony]]@queen@misc$yearOfBirth <- year
     }
   } else {
     stop("Argument x must be a Pop, Colony or MultiColony class object!")
diff --git a/R/Functions_L2_Colony.R b/R/Functions_L2_Colony.R
index 390e76f5..fcab0499 100644
--- a/R/Functions_L2_Colony.R
+++ b/R/Functions_L2_Colony.R
@@ -249,7 +249,7 @@ reQueen <- function(x, queen, removeVirginQueens = TRUE, simParamBee = NULL) {
 #' # nVirginQueens/nWorkers/nDrones will vary between function calls when a function is used
 #'
 #' # Queen's counters
-#' getMisc(getQueen(addWorkers(colony)))
+#' getQueen(addWorkers(colony))@misc
 #'
 #' # Add individuals to a MultiColony object
 #' apiary <- addWorkers(apiary)
@@ -472,7 +472,7 @@ addVirginQueens <- function(x, nInd = NULL, new = FALSE,
 #' nDrones(apiary)
 #'
 #' # Queen's counters
-#' getMisc(getQueen(buildUp(colony)))
+#' getQueen(buildUp(colony))@misc
 #' @export
 buildUp <- function(x, nWorkers = NULL, nDrones = NULL,
                     new = TRUE, exact = FALSE, resetEvents = FALSE,
diff --git a/man/addCastePop.Rd b/man/addCastePop.Rd
index 89e896e3..333abfc5 100644
--- a/man/addCastePop.Rd
+++ b/man/addCastePop.Rd
@@ -117,7 +117,7 @@ SP$nWorkers <- nWorkersPoisson
 # nVirginQueens/nWorkers/nDrones will vary between function calls when a function is used
 
 # Queen's counters
-getMisc(getQueen(addWorkers(colony)))
+getQueen(addWorkers(colony))@misc
 
 # Add individuals to a MultiColony object
 apiary <- addWorkers(apiary)
diff --git a/man/buildUp.Rd b/man/buildUp.Rd
index 5e280a04..4b6d03bc 100644
--- a/man/buildUp.Rd
+++ b/man/buildUp.Rd
@@ -118,5 +118,5 @@ nWorkers(apiary)
 nDrones(apiary)
 
 # Queen's counters
-getMisc(getQueen(buildUp(colony)))
+getQueen(buildUp(colony))@misc
 }
diff --git a/man/calcBeeGRMIbs.Rd b/man/calcBeeGRMIbs.Rd
index 2c32ebf4..f0fe324d 100644
--- a/man/calcBeeGRMIbs.Rd
+++ b/man/calcBeeGRMIbs.Rd
@@ -6,7 +6,7 @@
 \title{Calculate Genomic Relatedness Matrix (GRM) for honeybees from
   Identical By State genomic data}
 \usage{
-calcBeeGRMIbs(x, sex, alleleFreq = NULL)
+calcBeeGRMIbs(x, sex, alleleFreq = NULL, returnComponents = FALSE)
 
 calcBeeAlleleFreq(x, sex)
 }
@@ -20,10 +20,18 @@ missing values are allowed (this is not checked - you will get NAs!)}
 \code{x} - \code{"F"} for female and \code{"M"} for male}
 
 \item{alleleFreq}{numeric, vector of allele frequencies for the sites in
-\code{x}; if \code{NULL}, then \code{\link[SIMplyBee]{calcBeeAlleleFreq}} is used}
+\code{x}; if \code{NULL}, then \code{\link[SIMplyBee]{calcBeeAlleleFreq}}
+is used}
+
+\item{returnComponents}{logical, return GRM as well as the components used
+to compute it (useful for GWAS by GBLUP)}
 }
 \value{
-matrix of genomic relatedness coefficients
+if \code{returnComponents = FALSE} (default) return a matrix of
+  genomic relatedness coefficients; if \code{returnComponents = TRUE} return
+  a list with the GRM, centred genotype matrix, allele frequencies, and
+  scaling factor used to scale the crossproduct of centred genotype matrix
+  to get the GRM.
 }
 \description{
 Level 0 function that returns Genomic Relatedness Matrix (GRM)
@@ -56,12 +64,12 @@ sex <- getCasteSex(x = colony, collapse = TRUE)
 GRM <- calcBeeGRMIbs(x = geno, sex = sex)
 # You can visualise this matrix with the function image() from the package 'Matrix'
 
-#Look at the diagonal at the relationship matrix
+# Look at the diagonal at the relationship matrix
 x <- diag(GRM)
 hist(x)
 summary(x)
 
-#Look at the off-diagonal at the relationship matrix
+# Look at the off-diagonal at the relationship matrix
 x <- GRM[lower.tri(x = GRM, diag = FALSE)]
 hist(x)
 summary(x)
@@ -97,6 +105,10 @@ queenGeno <- getQueenSegSiteGeno(colony)
 calcBeeGRMIbs(x = rbind(queenGeno, pooledGenoW), sex = c("F","F"))
 # You can now compare how this compare to relationships between the queen
 # individual workers!
+
+# Return components
+calcBeeGRMIbs(x = rbind(queenGeno, pooledGenoW), sex = c("F","F"),
+              returnComponents = TRUE)
 }
 \references{
 Druet and Legarra (2020) Theoretical and empirical comparisons of
diff --git a/man/editCsdLocus.Rd b/man/editCsdLocus.Rd
index 080e039d..f713b776 100644
--- a/man/editCsdLocus.Rd
+++ b/man/editCsdLocus.Rd
@@ -32,3 +32,15 @@ Edits the csd locus in an entire population of individuals to
   recalculated to reflect the any changes due to editing, but other slots
   remain the same.
 }
+\examples{
+founderGenomes <- quickHaplo(nInd = 4, nChr = 1, segSites = 10)
+
+SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 2)
+tmp <- createVirginQueens(founderGenomes)
+getCsdAlleles(tmp)
+
+SP <- SimParamBee$new(founderGenomes, nCsdAlleles = 4)
+tmp <- createVirginQueens(founderGenomes)
+getCsdAlleles(tmp)
+
+}
diff --git a/man/getCasteSex.Rd b/man/getCasteSex.Rd
index db2f8c2a..e2bf1361 100644
--- a/man/getCasteSex.Rd
+++ b/man/getCasteSex.Rd
@@ -84,7 +84,5 @@ head(tmp)
 tail(tmp)
 }
 \seealso{
-\code{\link[SIMplyBee]{getCaste}}
-
 \code{\link[SIMplyBee]{getCastePop}} and \code{\link[SIMplyBee]{getCaste}}
 }
diff --git a/man/getIbdHaplo.Rd b/man/getIbdHaplo.Rd
index fd6dfe79..3386a215 100644
--- a/man/getIbdHaplo.Rd
+++ b/man/getIbdHaplo.Rd
@@ -165,5 +165,5 @@ getIbdHaplo(apiary, caste = "all")
 getIbdHaplo(apiary, caste = "all", collapse = TRUE)
 }
 \seealso{
-\code{\link[SIMplyBee]{getIbdHaplo}} and \code{\link[AlphaSimR]{pullIbdHaplo}}
+\code{\link[AlphaSimR]{pullIbdHaplo}}
 }
diff --git a/man/getMisc.Rd b/man/getMisc.Rd
deleted file mode 100644
index 60b342d8..00000000
--- a/man/getMisc.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/Functions_L0_auxilary.R
-\name{getMisc}
-\alias{getMisc}
-\title{Get miscellaneous information in a population}
-\usage{
-getMisc(x, node = NULL)
-}
-\arguments{
-\item{x}{\code{\link[AlphaSimR]{Pop-class}}}
-
-\item{node}{character, name of the node to get from the \code{x@misc} slot;
-if \code{NULL} the whole \code{x@misc} slot is returned}
-}
-\value{
-The \code{x@misc} slot or its nodes \code{x@misc[[*]][[node]]}
-}
-\description{
-Get miscellaneous information in a population
-}
diff --git a/man/getPooledGeno.Rd b/man/getPooledGeno.Rd
index c989d9e4..4cfcbf4f 100644
--- a/man/getPooledGeno.Rd
+++ b/man/getPooledGeno.Rd
@@ -4,21 +4,21 @@
 \alias{getPooledGeno}
 \title{Get a pooled genotype from true genotypes}
 \usage{
-getPooledGeno(x, type = NULL, sex = NULL)
+getPooledGeno(x, type = "mean", sex = NULL)
 }
 \arguments{
 \item{x}{matrix, true genotypes with individuals in rows and sites in columns}
 
-\item{type}{character, "mean" for average genotype or "count" for the counts
-of reference and alternative alleles}
+\item{type}{character, \code{"mean"} for average genotype (default) or
+\code{"count"} for the counts of reference and alternative alleles}
 
 \item{sex}{character, vector of "F" and "M" to denote the sex of individuals
 in \code{x}}
 }
 \value{
 a numeric vector with average allele dosage when \code{type = "mean"}
-  and a two-row matrix with the counts of reference (1st row) and
-  alternative (2nd row) alleles
+  or a two-row matrix with the counts of reference (1st row) and
+  alternative (2nd row) alleles when \code{type = "count"}
 }
 \description{
 Level 0 function that returns a pooled genotype from true
@@ -74,4 +74,9 @@ plot(
 
 # As an exercise you could repeat the above with different numbers of workers!
 
+# How to get pooled genotypes of workers across multiple colonies?
+tmp <- getSegSiteGeno(x = apiary, caste = "workers")
+(tmp2 = lapply(X = tmp, FUN = getPooledGeno)) # as a list of one row matrices
+t(sapply(X = tmp, FUN = getPooledGeno)) # as one matrix - option A
+do.call(what = rbind, args = tmp2) # as one matrix - option B
 }
diff --git a/man/getQtlGeno.Rd b/man/getQtlGeno.Rd
index 545a475b..16da4332 100644
--- a/man/getQtlGeno.Rd
+++ b/man/getQtlGeno.Rd
@@ -156,6 +156,8 @@ getQtlGeno(apiary, caste = "all")
 getQtlGeno(apiary, caste = "all", collapse = TRUE)
 }
 \seealso{
-\code{\link[SIMplyBee]{getQtlGeno}} and \code{\link[AlphaSimR]{pullQtlGeno}} as well as
+\code{\link[SIMplyBee]{getQtlHaplo}},
+  \code{\link[AlphaSimR]{pullQtlGeno}}, and
+  \code{\link[SIMplyBee]{getPooledGeno}}, as well as
   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 }
diff --git a/man/getQtlHaplo.Rd b/man/getQtlHaplo.Rd
index 353de61a..5cc4410f 100644
--- a/man/getQtlHaplo.Rd
+++ b/man/getQtlHaplo.Rd
@@ -173,6 +173,6 @@ getQtlHaplo(apiary, caste = "all", collapse = TRUE)
 
 }
 \seealso{
-\code{\link[SIMplyBee]{getQtlHaplo}} and \code{\link[AlphaSimR]{pullQtlHaplo}} as well as
+\code{\link[SIMplyBee]{getQtlGeno}} and \code{\link[AlphaSimR]{pullQtlHaplo}} as well as
   \code{vignette(topic = "QuantitativeGenetics", package = "SIMplyBee")}
 }
diff --git a/man/getSegSiteGeno.Rd b/man/getSegSiteGeno.Rd
index ed7300a4..d2d3c8db 100644
--- a/man/getSegSiteGeno.Rd
+++ b/man/getSegSiteGeno.Rd
@@ -148,5 +148,7 @@ getSegSiteGeno(apiary, caste = "all")
 getSegSiteGeno(apiary, caste = "all", collapse = TRUE)
 }
 \seealso{
-\code{\link[SIMplyBee]{getSegSiteGeno}} and \code{\link[AlphaSimR]{pullSegSiteGeno}}
+\code{\link[SIMplyBee]{getSegSiteHaplo}},
+  \code{\link[AlphaSimR]{pullSegSiteGeno}}, and
+  \code{\link[SIMplyBee]{getPooledGeno}}
 }
diff --git a/man/getSegSiteHaplo.Rd b/man/getSegSiteHaplo.Rd
index f24d27b5..e33b0919 100644
--- a/man/getSegSiteHaplo.Rd
+++ b/man/getSegSiteHaplo.Rd
@@ -163,5 +163,5 @@ getSegSiteHaplo(apiary, caste = "all")
 getSegSiteHaplo(apiary, caste = "all", collapse = TRUE)
 }
 \seealso{
-\code{\link[SIMplyBee]{getSegSiteHaplo}} and \code{\link[AlphaSimR]{pullSegSiteHaplo}}
+\code{\link[SIMplyBee]{getSegSiteGeno}} and \code{\link[AlphaSimR]{pullSegSiteHaplo}}
 }
diff --git a/man/getSnpGeno.Rd b/man/getSnpGeno.Rd
index ac68763e..7ba38324 100644
--- a/man/getSnpGeno.Rd
+++ b/man/getSnpGeno.Rd
@@ -161,5 +161,7 @@ getSnpGeno(apiary, caste = "all")
 getSnpGeno(apiary, caste = "all", collapse = TRUE)
 }
 \seealso{
-\code{\link[SIMplyBee]{getSnpGeno}} and \code{\link[AlphaSimR]{pullSnpGeno}}
+\code{\link[SIMplyBee]{getSnpHaplo}},
+ \code{\link[AlphaSimR]{pullSnpGeno}}, and
+ \code{\link[SIMplyBee]{getPooledGeno}}
 }
diff --git a/man/getSnpHaplo.Rd b/man/getSnpHaplo.Rd
index f424c5b8..0ebdcab3 100644
--- a/man/getSnpHaplo.Rd
+++ b/man/getSnpHaplo.Rd
@@ -172,5 +172,5 @@ getSnpHaplo(apiary, caste = "all", collapse = TRUE)
 
 }
 \seealso{
-\code{\link[SIMplyBee]{getSnpHaplo}} and \code{\link[AlphaSimR]{pullSnpHaplo}}
+\code{\link[SIMplyBee]{getSnpGeno}} and \code{\link[AlphaSimR]{pullSnpHaplo}}
 }
diff --git a/man/setMisc.Rd b/man/setMisc.Rd
deleted file mode 100644
index 0b2c0291..00000000
--- a/man/setMisc.Rd
+++ /dev/null
@@ -1,26 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/Functions_L0_auxilary.R
-\name{setMisc}
-\alias{setMisc}
-\title{Set miscellaneous information in a population}
-\usage{
-setMisc(x, node = NULL, value = NULL)
-}
-\arguments{
-\item{x}{\code{\link[AlphaSimR]{Pop-class}}}
-
-\item{node}{character, name of the node to set within the \code{x@misc} slot}
-
-\item{value, }{value to be saved into \code{x@misc[[*]][[node]]}; length of
-\code{value} should be equal to \code{nInd(x)}; if its length is 1, then
-it is repeated using \code{rep} (see examples)}
-}
-\value{
-\code{\link[AlphaSimR]{Pop-class}}
-}
-\description{
-Set miscellaneous information in a population
-}
-\details{
-A \code{NULL} in \code{value} is ignored
-}
diff --git a/vignettes/A_Honeybee_biology.Rmd b/vignettes/A_Honeybee_biology.Rmd
index bfc770d6..85cc1c3f 100644
--- a/vignettes/A_Honeybee_biology.Rmd
+++ b/vignettes/A_Honeybee_biology.Rmd
@@ -340,10 +340,10 @@ offspring that are expected to be homozygous based on queen's and father's *CSD*
 alleles. You can obtain `pHomBrood` and `nHomBrood` values with the
 corresponding `pHomBrood()` and `nHombrood()` functions that can be applied
 either on the queen (`Pop` class) or colony (`Colony` class) directly. You can
-obtain the entire `misc` slot with the `getMisc()` function.
+obtain the entire `misc` slot using:
 
 ```{r misc}
-getMisc(getQueen(colony))
+getQueen(colony)@misc
 ```
 
 Technically, in SIMplyBee we represent the *CSD* locus as a series of bi-allelic
@@ -412,7 +412,7 @@ population as well as the cumulative number of workers, drones, homozygous
 brood, and the expected proportion of homozygous brood.
 
 ```{r queens counters}
-getMisc(getQueen(inbredColony))
+getQueen(inbredColony)@misc
 ```
 
 # References