An extension of the ecologist’s guide to the animal model

For my current postdoc I have the opportunity to do some quantitative genetics with pedigree data in the animal model, something I’ve never done (having only used data from hemiclones/chromosome substitution line experiments before). I spent a couple of days learning how, and I turned to “An Ecologist’s Guide to the Animal Model” (Wilson et al 2009) for assistance. The tutorial included in the supplementary material (see revised version here, not the original supplementary material) proved to be an excellent and easy introduction. Being familiar with MCMCglmm, I used that as the environment for me to learn in, and I highly recommend this tutorial for anyone learning about the animal model.


For myself, I’m interested in the evolution of sex differences, which naturally calls for estimates of intersexual genetic correlations and G matrices. I decided to go beyond the tutorial and construct G matrices for the Gryphon data, which is included with the guide. To do this I had to generate sex-specific trait columns, and use a model with four response variables (female body weight [BWT_F], male body weight [BWT_M], female tarsus length [TAR_F], male tarsus length [TAR_M]). From this model I am able to estimate additive genetic variance, within-sex cross-trait covariance, cross-sex within-trait covariance, and cross-sex cross-trait covariance; these are the components of a G matrix. I can also derive correlations (e.g. intersexual genetic correlations, rMF), and within- and cross-sex heritability (see Lehtovaara et al 2013 Am Nat). The following annotated R script does all of this, and run for 1.3M iterations it produces the G matrix above (note that the numbers above the diagonal of G, the shaded cells, are the genetic correlations; the sex-specific additive genetic variances are in bold). I could also derive an M matrix; one with the variance and covariance explained by the maternal effects. Note, I haven’t worried about priors here (see the MCMCglmm documentation for advice), this was an exercise in extending the model and extracting intersexual genetic correlation /G matrices. It’s a work in progress – I would greatly appreciate feedback/corrections!

# Version 12-01-2017
# Load Package

# Load Data (download the data and change the filepaths to suit)
Data = read.table("\\gryphon.txt", header = T)
Ped = read.table("\\gryphonped.txt", header = T)

# Prep Data
colnames(Data)[1] = "animal"
Data$animal = as.factor(Data$animal)
Data$MOTHER = as.factor(Data$MOTHER)
Data$BYEAR = as.factor(Data$BYEAR)
Data$SEX = as.factor(Data$SEX)
Data$BWT = as.numeric(Data$BWT)
Data$TARSUS = as.numeric(Data$TARSUS)
for(x in 1:3) Ped[,x] = as.factor(Ped[,x])

# Make Sex-Specific Traits
Data$BWT_M = ifelse(Data$SEX == 2, Data$BWT, NA)
Data$BWT_F = ifelse(Data$SEX == 1, Data$BWT, NA) 
Data$TAR_M = ifelse(Data$SEX == 2, Data$TARSUS, NA)
Data$TAR_F = ifelse(Data$SEX == 1, Data$TARSUS, NA)

# Construct Priors (see Wilson et al 2009 and MCMCglmm documentation for guidance on priors)
phen.var = matrix(c(var(Data$BWT_F, na.rm = TRUE), 0, 0, 0,
 0, var(Data$BWT_M, na.rm = TRUE), 0, 0,
 0, 0, var(Data$TAR_F, na.rm = TRUE), 0,
 0, 0, 0, var(Data$TAR_M, na.rm = TRUE)
prior4.1= list(G = list( G1 = list(V = phen.var/4, n = 2),
 G2 = list(V = phen.var/4, n = 2),
 G3 = list(V = phen.var/4, n = 2)),
 R = list(V = phen.var/4, n = 2))

# Model
model4.1 = MCMCglmm(cbind(BWT_F, BWT_M, TAR_F, TAR_M) ~ trait-1,
 random = ~us(trait):animal + us(trait):BYEAR + us(trait):MOTHER,
 rcov = ~us(trait):units,
 family = rep("gaussian", times = 4),
 pedigree = Ped,
 data = Data,
prior = prior4.1)

# Additive Genetic Variance (Diagonals of G)
posterior.mode(VA_BWT_F <- model4.1$VCV[,"traitBWT_F:traitBWT_F.animal"])
posterior.mode(VA_BWT_M <- model4.1$VCV[,"traitBWT_M:traitBWT_M.animal"])
posterior.mode(VA_TAR_F <- model4.1$VCV[,"traitTAR_F:traitTAR_F.animal"])
posterior.mode(VA_TAR_M <- model4.1$VCV[,"traitTAR_M:traitTAR_M.animal"])

# Cross Trait Covariance (Off-diagonals of GM and GF)
posterior.mode(COV_BWT_F.TAR_F <- model4.1$VCV[,"traitBWT_F:traitTAR_F.animal"])
posterior.mode(COV_BWT_M.TAR_M <- model4.1$VCV[,"traitBWT_M:traitTAR_M.animal"])

# Cross Sex Covariance (Diagonals of B)
posterior.mode(COV_BWT_F.BWT_M <- model4.1$VCV[,"traitBWT_F:traitBWT_M.animal"])
posterior.mode(COV_TAR_F.TAR_M <- model4.1$VCV[,"traitTAR_F:traitTAR_M.animal"])

# Cross Trait Cross Sex Covariance (Off-diagonals of B)
posterior.mode(COV_BWT_F.TAR_M <- model4.1$VCV[,"traitBWT_F:traitTAR_M.animal"])
posterior.mode(COV_BWT_M.TAR_F <- model4.1$VCV[,"traitBWT_M:traitTAR_F.animal"])

# Cross Trait Genetic Correlation
posterior.mode(COR_BWT_F.TAR_F <- COV_BWT_F.TAR_F/sqrt(VA_BWT_F * VA_TAR_F))
posterior.mode(COR_BWT_M.TAR_M <- COV_BWT_M.TAR_M/sqrt(VA_BWT_M * VA_TAR_M))

# Cross Sex Genetic Correlation (AKA Intersexual Genetic Correlation, or rMF)
posterior.mode(COR_BWT_F.BWT_M <- COV_BWT_F.BWT_M/sqrt(VA_BWT_F * VA_BWT_M))
posterior.mode(COR_TAR_F.TAR_M <- COV_TAR_F.TAR_M/sqrt(VA_TAR_F * VA_TAR_M))

# Cross Trait Cross Sex Genetic Correlation
posterior.mode(COR_BWT_F.TAR_M <- COV_BWT_F.TAR_M/sqrt(VA_BWT_F * VA_TAR_M))
posterior.mode(COR_BWT_M.TAR_F <- COV_BWT_M.TAR_F/sqrt(VA_BWT_F * VA_TAR_M))

# Within-Sex Heritability (D or S = Dam or Sire, D or S = Daughter or Son)
posterior.mode(DDH2_BWT_F <- model4.1$VCV[,"traitBWT_F:traitBWT_F.animal"] / (model4.1$VCV[,"traitBWT_F:traitBWT_F.animal"] + model4.1$VCV[,"traitBWT_F:traitBWT_F.MOTHER"] + model4.1$VCV[,"traitBWT_F:traitBWT_F.BYEAR"] + model4.1$VCV[,"traitBWT_F:traitBWT_F.units"]))
posterior.mode(SSH2_BWT_M <- model4.1$VCV[,"traitBWT_M:traitBWT_M.animal"] / (model4.1$VCV[,"traitBWT_M:traitBWT_M.animal"] + model4.1$VCV[,"traitBWT_M:traitBWT_M.MOTHER"] + model4.1$VCV[,"traitBWT_M:traitBWT_M.BYEAR"] + model4.1$VCV[,"traitBWT_M:traitBWT_M.units"]))
posterior.mode(DDH2_TAR_F <- model4.1$VCV[,"traitTAR_F:traitTAR_F.animal"] / (model4.1$VCV[,"traitTAR_F:traitTAR_F.animal"] + model4.1$VCV[,"traitTAR_F:traitTAR_F.MOTHER"] + model4.1$VCV[,"traitTAR_F:traitTAR_F.BYEAR"] + model4.1$VCV[,"traitTAR_F:traitTAR_F.units"]))
posterior.mode(SSH2_TAR_M <- model4.1$VCV[,"traitTAR_M:traitTAR_M.animal"] / (model4.1$VCV[,"traitTAR_M:traitTAR_M.animal"] + model4.1$VCV[,"traitTAR_M:traitTAR_M.MOTHER"] + model4.1$VCV[,"traitTAR_M:traitTAR_M.BYEAR"] + model4.1$VCV[,"traitTAR_M:traitTAR_M.units"]))

# Cross-Sex Heritability (D or S = Dam or Sire, D or S = Daughter or Son)
posterior.mode(DSH2_BWT <- DDH2_BWT_F * COR_BWT_F.BWT_M)
posterior.mode(SDH2_BWT <- SSH2_BWT_M * COR_BWT_F.BWT_M)
posterior.mode(DSH2_TAR <- DDH2_TAR_F * COR_TAR_F.TAR_M)
posterior.mode(SDH2_TAR <- SSH2_TAR_M * COR_TAR_F.TAR_M)

# G matrix
G = matrix(c( posterior.mode(VA_BWT_F), posterior.mode(COR_BWT_F.TAR_F), posterior.mode(COR_BWT_F.BWT_M), posterior.mode(COR_BWT_M.TAR_F),
 posterior.mode(COV_BWT_F.TAR_F), posterior.mode(VA_TAR_F), posterior.mode(COR_BWT_F.TAR_M), posterior.mode(COR_TAR_F.TAR_M),
 posterior.mode(COV_BWT_F.BWT_M), posterior.mode(COV_BWT_M.TAR_F), posterior.mode(VA_BWT_M), posterior.mode(COR_BWT_M.TAR_M),
 posterior.mode(COV_BWT_F.TAR_M), posterior.mode(COV_TAR_F.TAR_M), posterior.mode(COV_BWT_M.TAR_M), posterior.mode(VA_TAR_M))
 , ncol = 4, byrow=T)

colnames(G) = rownames(G) = c("BWT_F", "TAR_F", "BWT_M", "TAR_M")


New Paper: Constrained Evolution

Directional selection should erode genetic variance, but covariance among traits and antagonistic selection in different contexts could maintain it. In a recent project, led by Clarissa House and linked back to my MSc work, we used a quantitative genetic approach to look at the evolution of sex comb morphology in fruit flies, studying selection in pre- and post-copulatory contexts. While we show genetic variance within the traits and genetic correlations among the traits, it appeared that a lack of directional selection on the multivariate phenotype is in fact placing limits on adaptive evolution, allowing genetic variation in sex comb morphology to persist. The paper is now available via the Journal of Evolutionary Biology.




ESEB 2017 Symposium

The ESEB 2017 Scientific Committee have confirmed that Fiona Ingleby and I will be hosting a symposium at the next meeting in Groningen. The symposium, Modern quantitative genetics and the study of adaptation, will give a broad perspective on recent advances in quantitative genetic approaches (for example: the use of omics, hemiclone and chromosome substitution lines, methods to estimate large G-matrices), and how these can be used to study a number of evolutionary processes with adaptation at their core (for example: the evolution of sex differences, host-parasite/predator-prey co-evolutionary dynamics, evolutionary medicine, conservation strategy and adaptation to climate change).


We are very excited that Mark Blows and Anne Charmantier will be joining us as our plenary speakers and offering their considerable expertise on multivariate evolution, quantitative genetics in natural populations, and adaptation. We aim to welcome abstracts from researchers using modern quantitative genetic methods to tackle a variety of problems at the fore of modern evolutionary biology, with a strong emphasis on encouraging a cross-disciplinary exchange of knowledge and ideas.

Looking forward to ESEB2017!

New Paper: X-linkage of Sex-Specific Genetic Variance

The third paper from my PhD has been published in G3: Genes, Genomes, Genetics. In this paper we developed a new method to separately estimate sex-specific genetic variance and cross-sex genetic covariance of the X chromosomes and major autosomes in Drosophila melanogaster, and assayed lifespan and ageing.


The population genetics of the X chromosome differ to those of the autosomes because they occur as a single copy in males, just like in humans. This means that sexually antagonistic and, consequently, sex-specific mutations may more readily accumulate on the X, a prediction supported by our findings. We show that the cross-sex genetic correlation is lower for lifespan on the X chromosome, suggesting a potential for the X to play a special role in the evolution of sex differences.


Genetic Drift Simulator For R

Here is a program to simulate genetic drift in R (code pasted below), specifically one capable of handling unequal sex ratio. Sex ratio directly affects the effective population size, and this simulation can be used to show that it is the effective population size (NE), rather than the census population size (NA), that determines the population genetic behaviour of a population.

Here are some demonstration results, where 10 replicate populations are simulated with NE = 360, but very different actual population sizes, and an initial allele frequency of 0.5. The left hand simulations have unequal sex ratios, while the sex ratio is equal in the right hand simulations. Both sets of simulations behave similarly, showing that effective population size, rather than actual population size, is the determinant of population genetic behaviour.


#####~ Genetic Drift Simulator ~#####
#######~  Robert M. Griffin  ~#######
#####~ Updated March 22nd 2016 ~#####
# Pop   = Replicate populations
# Gen   = Generations
# NM    = Male population size
# NF    = Female population size
# P     = Frequency of focal allele
# histo = Optional histogram
GenDriftSim = function(Pop = Pop, Gen = Gen, NM, NF, P, histo = "yes"){
 P = (2*(NM+NF))*P
 NE = round((4*NM*NF)/(NM+NF),0)
 SR = round(NM/NF,2)
 Na = NM+NF
 plot(c(0,0),type = "n", main = bquote('N'[M]~'/ N'[F]~'='~.(SR)*', N'[A]~'='~.(Na)*', N'[E]~'='~.(NE)), cex.main = 1, 
 xlim = c(1,Gen), ylim=c(0,1), xlab = "Generations", ylab = "Fequency of focal allele")
 for (i in 1:Pop){
 N = NM+NF
 startA = as.vector(c(rep(1, times = P),rep(0, times = (2*N)-P)))
 Population = matrix(c(
 c(sample(startA, size = 2*N, replace = FALSE)),
 c(rep("M", times = NM), rep("F", times = NF))),
 ncol = 3)
 SimResults[(Gen*i)+1-Gen, 3] <<- sum(as.numeric(Population[,1:2]))/(N*2)
 for(j in 1:(Gen-1)){
 Population = matrix(c(
 c(sample(sample(Population[(1:NM),1:2], replace = TRUE),N, replace = TRUE)),
 c(sample(sample(Population[(1+NM):N,1:2], replace = TRUE),N, replace = TRUE)),
 c(rep("M", times = NM), rep("F", times = NF))), ncol = 3)
 SimResults[(Gen*i)+1+j-Gen, 3] <<- sum(as.numeric(Population[,1:2]))/(N*2)
 s = (i*Gen)-Gen+1; e = i*Gen
 r = as.vector(SimResults[s:e, 3])
 points(r~c(1:Gen), type = "l")
 if(histo == "yes"){SimResults[,1] = rep(1:Pop, each = Gen)
 SimResults[,2] = rep(1:Gen, times = Pop)
 hist(SimResults[,3][SimResults[,2]==Gen], xlim = c(0,1), cex.main = 1, main = bquote('N'[M]~'/ N'[F]~'='~.(SR)*', N'[A]~'='~.(Na)*', N'[E]~'='~.(NE)), xlab = paste0("Frequency of focal allele after ",Gen," Generations"))

Pop = 10
Gen = 25
P = 0.5
histo = "no"

SimResults = matrix(data = NA, ncol = 3, nrow = Gen*Pop)
GenDriftSim(Pop = Pop, Gen = Gen, NM = 100, NF = 900, P = P, histo = histo)
GenDriftSim(Pop = Pop, Gen = Gen, NM = 180, NF = 180, P = P, histo = "histo")

Postdoctoral Position

I am now a postdoctoral researcher on natural selection in contemporary human populations at Turku University in Virpi Lummaa’s research group! Virpi’s group studies life history evolution in Humans and Elephants using demographic data. I am really excited to be joining the group as they move from Sheffield to Turku, and look forward to trying my hand at quantitative genetics with pedigree data! News on projects will come soon.


Göteborg University Talk

I gave a presentation at Göteborg University Department of Zoology today, titled The Evolution of Sexual Dimorphism: Quantitative Genetic Explorations of the Drosophila Genome, which was based on my PhD work at Uppsala University. It was really good to present the complex elements of my work to a diverse research department, explaining the concept of adaptation and constraint in multivariate space over just 45 minutes is quite a challenge!


Thanks to Lotta Kvarnemo and the Department of Zoology for inviting me to speak today and the good discussion that followed.