#--------------------------------------------------------------# # Data analysis for the manuscript: # "The Potential for Genotype-by-Environment Interactions to # Maintain Genetic Variation in a Model Legume–Rhizobia Mutualism" # Plant Communications (Volume 1, Issue 6) # 9 November 2020 # Priya Vaidya and John Stinchcombe # Corresponding author: Priya Vaidya (p.vaidya@mail.utoronto.ca) #--------------------------------------------------------------# #--------------------------------------------------------------# ##### Prepare workspace ##### #--------------------------------------------------------------# # Load libraries library(lme4) library(car) library(Hmisc) library(emmeans) library(corrplot) library(lattice) library(ggplot2) library(dplyr) library(reshape) # Load script with custom functions to check assumptions of linear models (see Wood et al., 2018) source("/Users/priyavaidya/Desktop/Dryad/checkAssumptions.R") # Set contrasts (sum-to-zero rather than R's default treatment contrasts to account for unequal sample sizes across treatments) # http://rstudio-pubs-static.s3.amazonaws.com/65059_586f394d8eb84f84b1baaf56ffb6b47f.html options(contrasts=c("contr.sum", "contr.poly")) # Read in data data <- read.csv("/Users/priyavaidya/Desktop/Dryad/Raw Data - Vaidya & Stinchcombe Plant Comm.csv") # Remove Controls from dataset data = subset(data, Short.ID!="Control") # Create unique identifier that includes Treatment and Genotype of an individual plant data$TG = paste(data$Treatment, data$Genotype, sep="_") View(data) #--------------------------------------------------------------# #--------------------------------------------------------------# # DATA ANALYSIS #--------------------------------------------------------------# #--------------------------------------------------------------# ##### Effects of Light Treatment, Plant Genotype, and their interaction (GxE) on log-transformed Plant Performance and Mutualism-Related Traits ##### #--------------------------------------------------------------# # Leaf number # Model 1: Plant genotype and GxE included as random variables leaf1 <- lmer(log(leaves+1)~Treatment + (1|Genotype) + (1|Treatment:Genotype) + (1|Block) + (1|Block:Treatment) + (1|Block:Genotype), data, control=lmerControl(optimizer="bobyqa")) check.assumptions(leaf1) # log transformed data meet assumptions Anova(leaf1, type=3, test="F") summary(leaf1) # Likelihood ratio tests for effect of GxE (Model 2) and Plant genotype (Model 3) on nodule number # Model 2: GxE excluded from Model 1 leaf2 = update(leaf1, .~. - (1|Treatment:Genotype)) # Model 3: Plant genotype excluded from Model 1 leaf3 = update(leaf1, .~. - (1|Genotype)) # NOTE: one-sided test = **p-values must be halved** anova(leaf1, leaf2) anova(leaf1, leaf3) # Nodule number # Include total biomass as covariate to account for whole plant response to light availability # Model 1: Plant genotype and GxE included as random variables nod1 <- lmer(log(nodules+1)~Treatment + total.biomass + (1|Genotype) + (1|Treatment:Genotype) + (1|Block) + (1|Block:Treatment) + (1|Block:Genotype), data, control=lmerControl(optimizer="bobyqa")) check.assumptions(nod1) # Assumptions met summary(nod1) Anova(nod1, type=3, test="F") # Likelihood ratio tests for effect of GxE (Model 2) and Plant genotype (Model 3) on nodule number # Model 2: GxE excluded from Model 1 nod2 = update(nod1, .~. - (1|Treatment:Genotype)) # Model 3: Plant genotype excluded from Model 1 nod3 = update(nod1, .~. - (1|Genotype)) # NOTE: one-sided test = **p-values must be halved** anova(nod1, nod2) anova(nod1, nod3) # Aboveground (shoot) biomass # Model 1: Plant genotype and GxE included as random variables s1 <- lmer(log(shoot.biomass+1)~Treatment + (1|Genotype) + (1|Treatment:Genotype) + (1|Block) + (1|Block:Treatment) + (1|Block:Genotype), data, control=lmerControl(optimizer="bobyqa")) check.assumptions(s1) #Assumptions met Anova(s1, type=3, test="F") summary(s1) # Likelihood ratio tests for effect of GxE (Model 2) and Plant genotype (Model 3) on nodule number # Model 2: GxE excluded from Model 1 s2=update(s1, .~. - (1|Treatment:Genotype)) # Model 3: Plant genotype excluded from Model 1 s3=update(s1, .~. - (1|Genotype)) # NOTE: one-sided test = **p-values must be halved** anova(s1, s2) anova(s1,s3) # Belowground (root) biomass # Model 1: Plant genotype and GxE included as random variables r1 <- lmer(log(root.biomass+1)~Treatment + (1|Genotype) + (1|Treatment:Genotype) + (1|Block) + (1|Block:Treatment) + (1|Block:Genotype), data, control=lmerControl(optimizer="bobyqa")) check.assumptions(r1) #Assumptions met Anova(r1, type=3, test="F") summary(r1) # Likelihood ratio tests for effect of GxE (Model 2) and Plant genotype (Model 3) on nodule number # Model 2: GxE excluded from Model 1 r2 <- update(r1, .~. - (1|Treatment:Genotype)) # Model 3: Plant genotype excluded from Model 1 r3 <- update(r1, .~. - (1|Genotype)) # NOTE: one-sided test = **p-values must be halved** anova(r1,r2) anova(r1,r3) # Aboveground:Belowground (shoot:root) biomass # Model 1: Plant genotype and GxE included as random variables sr1 <- lmer(log(shoot.root+1)~Treatment + (1|Genotype) + (1|Treatment:Genotype) + (1|Block) + (1|Block:Treatment) + (1|Block:Genotype), data, control=lmerControl(optimizer="bobyqa")) check.assumptions(sr1) Anova(sr1, type=3, test="F") summart(sr1) # Likelihood ratio tests for effect of GxE (Model 2) and Plant genotype (Model 3) on nodule number # Model 2: GxE excluded from Model 1 sr2 <- update(sr1, .~. - (1|Treatment:Genotype)) # Model 3: Plant genotype excluded from Model 1 sr3 <- update(sr1, .~. - (1|Genotype)) # NOTE: one-sided test = **p-values must be halved** anova(sr1,sr2) anova(sr1,sr3) #--------------------------------------------------------------# ##### Genetic Correlations Among Traits and Across Treatments ##### #--------------------------------------------------------------# # Create line means of raw data for each plant performance and mutualism-related trait #Leaf number leafrawmeans <- data.frame(group_by(data, TG, Treatment, Genotype) %>% dplyr::summarize(rawmean = mean(leaves))) leafrawmeans$Model = "leafrawmeans" #Nodule number nodrawmeans <- data.frame(group_by(data, TG, Treatment, Genotype) %>% dplyr::summarize(rawmean = mean(nodules))) nodrawmeans$Model = "nodrawmeans" #Aboveground (shoot) biomass srawmeans <- data.frame(group_by(data, TG, Treatment, Genotype) %>% dplyr::summarize(rawmean = mean(shoot.biomass))) srawmeans$Model = "srawmeans" #Belowground (root) biomass rrawmeans <- data.frame(group_by(data, TG, Treatment, Genotype) %>% dplyr::summarize(rawmean = mean(root.biomass))) rrawmeans$Model = "rrawmeans" # Aboveground:Belowground (shoot:root) biomass srrawmeans <- data.frame(group_by(data, TG, Treatment, Genotype) %>% dplyr::summarize(rawmean = mean(shoot.root))) srrawmeans$Model = "srrawmeans" # Combine raw line means of each trait into one data frame rawmean <- rbind(leafrawmeans, nodrawmeans, srawmeans, rrawmeans,srrawmeans) rawmean <- reshape(rawmean, direction= "wide", timevar="Model", v.names="rawmean", idvar="TG") names(rawmean)[4:8] <- c("Leaf #", "Nodule #", "Aboveground Biomass", "Belowground Biomass", "Aboveground:Belowground Ratio") # Use raw line means to create correlation plot #Separate raw line means data frame into Ambient and Shade line means for each trait treat <- reshape(rawmean, direction= "wide", timevar="Treatment", idvar="Genotype", drop=c("TG")) names(treat)[2:11] <- c("Ambient Leaf #", "Ambient Nodule #", "Ambient Aboveground Biomass", "Ambient Belowground Biomass", "Ambient Aboveground:Belowground Ratio","Shade Leaf #", "Shade Nodule #", "Shade Aboveground Biomass", "Shade Belowground Biomass", "Shade Aboveground:Belowground Ratio") treat <- treat[-(1)] # Colour coding mycolors <- rep(NA,10) names(mycolors) <- names(treat) mycolors[1:5] <- '#e2b07a' #Ambient colour mycolors[6:10] <- '#66a8a1' #Shade colour #Calculate Pearson correlation coefficients (r) and p-values treatcorr <- rcorr(as.matrix(treat),type="pearson") #Graph the correlation matrix #NOTE: The orientation of the correlation matrix is inverted in the original article, # but the specific correlation coefficients and p-values remain consistent corrplot(treatcorr$r,type="upper",tl.col = mycolors,method="ellipse", tl.srt=45, tl.cex = 0.6, addCoef.col="black", cl.offset = 0.3, cl.align.text = "l", p.mat=treatcorr$P, sig.level=0.05, insig = "blank", number.cex = 0.85) #--------------------------------------------------------------# ##### Significant GxE interactions due to changes in rank order versus changes in the magnitude of genetic variance ##### #--------------------------------------------------------------# #NOTE: Only traits that showed a significant (p<0.05) GxE interaction were assessed for Rank vs. Variance #NOTE: Calculations call for square-root of the genotype variance; use genotype SD #Nodule number # Genetic variance (Vg) in Shade treatment nodshade <- lmer(log(nodules+1)~ total.biomass + (1|Genotype) + (1|Block), subset(data, Treatment!="Ambient"), control=lmerControl(optimizer="bobyqa")) summary(nodshade) # Genetic variance (Vg) in Ambient treatment nodamb <- lmer(log(nodules+1)~ total.biomass + (1|Genotype) + (1|Block), subset(data, Treatment!="Shade"), control=lmerControl(optimizer="bobyqa")) summary(nodamb) # Pearson correlation coefficient (r) between Shade and Ambient treatments using log transformed line means nodmeans <- data.frame(group_by(data, TG, Treatment, Genotype) %>% dplyr::summarize(rawmean = mean(log(nodules+1)))) nodmeans <- reshape(nodmeans, direction= "wide", timevar="Treatment", idvar="Genotype", drop=c("TG")) cor.test(nodmeans$rawmean.Ambient,nodmeans$rawmean.Shade, method="pearson") # Change in rank (Vrank) = 2*Shade_Vg*Ambient_Vg*(1-r) 2*0.1772*0.2528*(1-0.6821529) # = 0.02847666 # Change in variance (Vvariance) = (Shade_Vg - Ambient_Vg)^2 (0.1772 - 0.2528)^2 # = 0.00571536 # Vgxe due to change in rank (Vrank/(Vrank + Vvariance)) 0.02847666/(0.02847666+0.00571536) # = 0.8328452 or 83.3% #Aboveground (shoot) biomass # Genetic variance (Vg) in Shade treatment sshade <- lmer(log(shoot.biomass+1)~ (1|Genotype) + (1|Block), subset(data,Treatment!="Ambient"), control=lmerControl(optimizer="bobyqa")) summary(sshade) # Genetic variance (Vg) in Ambient treatment samb <- lmer(log(shoot.biomass+1)~ (1|Genotype) + (1|Block), subset(data,Treatment!="Shade"), control=lmerControl(optimizer="bobyqa")) summary(samb) # Pearson correlation coefficient (r) between Shade and Ambient treatments using log transformed line means shootmeans <- data.frame(group_by(data, TG, Treatment, Genotype) %>% dplyr::summarize(rawmean = mean(log(shoot.biomass+1)))) shootmeans <- reshape(shootmeans, direction= "wide", timevar="Treatment", idvar="Genotype", drop=c("TG")) cor.test(shootmeans$rawmean.Ambient,shootmeans$rawmean.Shade, method="pearson") # Change in rank (Vrank) = 2*Shade_Vg*Ambient_Vg*(1-r) 2*0.02616*0.09470*(1-0.5615821) # = 0.002172231 # Change in variance (Vvariance) = (Shade_Vg - Ambient_Vg)^2 (0.02616 - 0.09470)^2 # = 0.004697732 # Vgxe due to change in rank (Vrank/(Vrank + Vvariance)) 0.002172231/(0.002172231+0.004697732) # = 0.3161925 or 31.6% #Belowground (root) biomass # Genetic variance (Vg) in Shade treatment rshade <- lmer(log(root.biomass+1)~ (1|Genotype) + (1|Block), subset(data,Treatment!="Ambient"), control=lmerControl(optimizer="bobyqa")) summary(rshade) # Genetic variance (Vg) in Ambient treatment ramb <- lmer(log(root.biomass+1)~ (1|Genotype) + (1|Block), subset(data,Treatment!="Shade"), control=lmerControl(optimizer="bobyqa")) summary(ramb) # Pearson correlation coefficient (r) between Shade and Ambient treatments using log transformed line means rootmeans <- data.frame(group_by(data, TG, Treatment, Genotype) %>% dplyr::summarize(rawmean = mean(log(root.biomass+1)))) rootmeans <- reshape(rootmeans, direction= "wide", timevar="Treatment", idvar="Genotype", drop=c("TG")) cor.test(rootmeans$rawmean.Ambient,rootmeans$rawmean.Shade, method="pearson") # Change in rank (Vrank) = 2*Shade_Vg*Ambient_Vg*(1-r) 2*0.009006*0.04402*(1-0.779015) # = 0.0001752164 # Change in variance (Vvariance) = (Shade_Vg - Ambient_Vg)^2 (0.009006 - 0.04402)^2 # = 0.00122598 # Vgxe due to change in rank (Vrank/(Vrank + Vvariance)) 0.0001752164/(0.0001752164+0.00122598) # = 0.1250477 or 12.5% #--------------------------------------------------------------# # SELECTION ANALYSIS #--------------------------------------------------------------# # Obtain raw line means of the response variable and focal traits (excluding zeros) # Aboveground (shoot) biomass - response variable (fitness proxy) slinemean <- data.frame(group_by(subset(data, shoot.biomass>0),TG, Treatment, Genotype) %>% dplyr::summarize(linemean = mean(shoot.biomass))) slinemean$Model = "slinemean" # Nodule number (focal mutualistic trait) nodlinemean <- data.frame(group_by(subset(data, nodules>0),TG, Treatment, Genotype) %>% dplyr::summarize(linemean = mean(nodules))) nodlinemean$Model = "nodlinemean" # Leaf number (plant performance trait) leaflinemean <- data.frame(group_by(subset(data, leaves>0), TG, Treatment, Genotype) %>% dplyr::summarize(linemean = mean(leaves))) leaflinemean$Model = "leaflinemean" # Belowground (root) biomass (plant performance trait) rlinemean <- data.frame(group_by(subset(data, root.biomass>0),TG, Treatment, Genotype) %>% dplyr::summarize(linemean = mean(root.biomass))) rlinemean$Model = "rlinemean" # Combine trait line means into one data frame linemeans <- rbind(leaflinemean, nodlinemean, slinemean, rlinemean) linemeans <- reshape(linemeans, direction= "wide", timevar="Model", v.names="linemean", idvar="TG") names(linemeans)[4:7] <- c("leaflinemean", "nodlinemean", "slinemean", "rlinemean") #--------------------------------------------------------------# ##### Absolute fitness ##### #--------------------------------------------------------------# # No relativization of fitness proxy or standardization of focal traits # Assessment of relationship between fitness and focal traits as measured # Selection-by-treatment (SxE) analysis absfit <- lm(slinemean ~ nodlinemean*Treatment + leaflinemean*Treatment + rlinemean*Treatment, linemeans) Anova(absfit,type=3,test="F") # Estimates of selection gradient (beta) within each light treatment absfitamb <- lm(slinemean ~ nodlinemean + leaflinemean + rlinemean, subset(linemeans,Treatment!="Shade")) summary(absfitamb) absfitshade <- lm(slinemean ~ nodlinemean + leaflinemean + rlinemean, subset(linemeans,Treatment!="Ambient")) summary(absfitshade) #--------------------------------------------------------------# ##### Global scale (relativizing fitness and standardizing traits across treatments) ###### #--------------------------------------------------------------# # Calculate relative fitness (aboveground biomass) using the mean from across both light treatments relshoot <- linemeans$slinemean/mean(linemeans$slinemean) # Add relative fitness to data frame sel <- cbind(linemeans,relshoot) # Calculate standardized focal traits using means and standard deviations from across both light treatments sel <- cbind(sel[,-c(4:7)], rescaler(sel[,c(4,5,7)],"sd")) names(sel)[5:7] <- c("sdleaf", "sdnod", "sdroot") # Selection-by-treatment (SxE) analysis global <- lm(relshoot ~ Treatment*sdnod+ sdleaf*Treatment + sdroot*Treatment, sel) Anova(global, type=3, test="F") # Estimates of selection gradient (beta) within each light treatment globamb <- lm(relshoot ~ sdnod + sdleaf + sdroot, data = subset(sel, Treatment!="Shade")) summary(globamb) globshade <- lm(relshoot ~ sdnod + sdleaf + sdroot, data = subset(sel, Treatment!="Ambient")) summary(globshade) #--------------------------------------------------------------# ##### Local scale (relativizing fitness and standardizing traits within each treatments) ###### #--------------------------------------------------------------# # Calculate relative fitness (aboveground biomass) and standardized focal traits # using the means from within each light treatment #AMBIENT # Data frame of line means from the Ambient treatment selamb <- subset(linemeans, Treatment!="Shade") # Relative shoot biomass in Ambient ambrelshoot <- selamb$slinemean/mean(selamb$slinemean,na.rm=TRUE) # Add relative fitness to Ambient data frame selamb <- cbind(selamb,ambrelshoot) # Calculate standardized focal traits using means and standard deviations in the Ambient treatment selamb <- cbind(selamb[,-c(4:7)], rescaler(selamb[,c(4,5,7)],"sd")) # Standardized traits in Ambient names(selamb)[4:7] <- c("relshoot", "sdleaf", "sdnod", "sdroot") # SHADE # Data frame of line means from the Shade treatment selshade <- subset(linemeans, Treatment!="Ambient") # Relative shoot biomass in Shade shaderelshoot <- selshade$slinemean/mean(selshade$slinemean,na.rm=TRUE) # Add relative fitness to Shade data frame selshade <- cbind(selshade,shaderelshoot) # Calculate standardized focal traits using means and standard deviations in the Shade treatment selshade <- cbind(selshade[,-c(4:7)], rescaler(selshade[,c(4,5,7)],"sd")) names(selshade)[4:7] <- c("relshoot","sdleaf", "sdnod", "sdroot") View(selshade) # Combine Shade and Ambient data frames into one data frame sel2 <- rbind(selshade,selamb) # Selection-by-treatment (SxE) analysis selgrad2 <- lm(relshoot ~ sdnod*Treatment + sdleaf*Treatment + sdroot*Treatment, sel2) Anova(selgrad2, type=3, test="F") # Estimates of selection gradient (beta) within each light treatment sel2amb <- lm(relshoot ~ sdnod + sdleaf + sdroot, selamb) summary(sel2amb) sel2shade <- lm(relshoot ~ sdnod + sdleaf + sdroot, selshade) summary(sel2shade)