#### FORAGING BEHAVIOR AND EXTENDED PHENOTYPE INDEPENDENTLY AFFECT FORAGING SUCCESS IN BLACK WIDOW SPIDERS #### #### DIRIENZO N, SCHRAFT HA, MONTIGLIO PO, BRADLEY CT, DORNHAUS A #### ### Setup ### # remove all objects from workspace so we can start with clean slate rm(list = ls()) # load required packages library(lme4) # fitting mixed models library(dplyr) # data manipulation library(rptR) # calculate repeatabilities library(MuMIn) ## model comparison, calculating R2's library(ggplot2) # for plotting ### Data wrangling ### ## Pre data ### # This is the dataframe with the initial web assay and behavioral assay (i.e. before mesocosm trials). # load data behav <- read.csv("behav.csv", na = 'NA') #convert distance to numeric behav$distance.n <- ifelse(behav$distance == "f", 3, ifelse(behav$distance == "m", 2, 1)) # convert to tidy frame behav.df <- tbl_df(behav) # group by ID behav.ID <- behav.df %>% group_by(ID, web) # sum up attack number per group attack.number <- summarize(behav.ID, attack.number = sum(attack), retreat.perc = sum(retreat, na.rm=TRUE)/sum(attack)) # merge back into behav behav.df2 <- inner_join(behav.df, attack.number, by = c("ID", "web")) # center weight behav.df2 <- behav.df2 %>% mutate(weight.c = (weight - mean(weight))/sd(weight)) # make a DF of just the first observations so we only have one per web web.df <- filter(behav.df2, day == 1, obs.day == 1) # make dataframe to have rprR try to calculate repeatability of attack.number web.df2 <- web.df web.df2$attack.number <- attack.number$attack.number ## Mesocosm data ## # This is the data from the mesocosm trials. # import mesocosm data meso <- read.csv("meso.csv", na = 'NA') # to tibble meso.df <- tbl_df(meso) # group by ID meso.df <- meso %>% group_by(ID, web) # merge in attack number meso.df2 <- inner_join(meso.df, attack.number, by = c("ID", "web")) # center/scale meso.df2$weight.c <- scale(meso.df2$weight) meso.df2$web.ID.weight.c <- scale(meso.df2$web.ID.weight) # make time a factor meso.df2$time.factor <- as.factor(meso.df2$time) # calculate difference in gumfooted lines meso.df2$gumfooted.difference <- meso.df2$gumfooted.own - meso.df2$gumfooted # add OLRE meso.df2$INDV <- 1:nrow(meso.df2) meso.df2$gumfooted.c <- scale(meso.df2$gumfooted) meso.df2$gumfooted.own.c <- scale(meso.df2$gumfooted.own) meso.df2$time.c <- scale(meso.df2$time) meso.df2$attack.number.c <- scale(meso.df2$attack.number) meso.df2$structural.c <- scale(meso.df2$structural) meso.df2$structural.own.c <- scale(meso.df2$structural.own) meso.df2$gumfooted.diff <- meso.df2$gumfooted - meso.df2$gumfooted.own meso.df2$gumfooted.diff.c <- scale(meso.df2$gumfooted.diff) meso.df2$structural.diff <- meso.df2$structural - meso.df2$structural.own meso.df2$structural.diff.c <- scale(meso.df2$structural.diff) meso.df2$webID <- meso.df2$web.ID meso.df2$prey.number <- 5 meso.df2$prey.por <- meso.df2$interval.ants.eaten/5 # meso.df2$attack.half <- ifelse(meso.df2$attack.number >= 6.3, 1, 0) meso.all <- meso.df2 meso.swap <- meso.df2[meso.df2$swapped == 1,] meso.own <- meso.df2[meso.df2$swapped == 0,] # get data ready for repeatability / phenotypic correlation calculations meso.all.1 <- meso.all[meso.all$time == 1, ] meso.all.1$olre <- seq(1:nrow(meso.all.1)) # add OLRE for global model temporal autocorrelation meso.all$olre <- seq(1:nrow(meso.all)) # add OLRE for swap model temporal autocorrelation meso.swap$olre <- seq(1:nrow(meso.swap)) ## Post data ## # This is the web structure data taken after the mesocosm trial. It shows how much webs changed during the mesocosm trial. # load data post <- read.csv("post.csv", na = 'NA') # to tibble post.df <- tbl_df(post) post.df2 <- inner_join(post.df, attack.number, by = c("ID", "web")) # extract total prey captured from 96 hour interval in meso.df2 # filter just the last day and then select just the columns we care about prey.captured <- filter(meso.df2, time == 96 ) prey.captured <- select_(prey.captured, "ID", "web", "total.ants.eaten", "total.crickets.eaten", "attack.number.c", "gumfooted.c", "structural.c", "gumfooted.diff.c", "structural.diff.c", "weight.c", "gumfooted.diff", "structural.diff") # merge in prey capture data to our post mesocosm data frame post.df3 <- inner_join(post.df2, prey.captured, by = c("ID", "web")) # Change in web structure over time post.df3$structural.change <- post.df3$structural.after - post.df3$structural post.df3$gumfooted.change <- post.df3$gumfooted.after - post.df3$gumfooted # difference in web from own and how it affects prey capture post.df3$structural.from.own <- post.df3$structural - post.df3$structural.own post.df3$gumfooted.from.own <- post.df3$gumfooted - post.df3$gumfooted.own # also add changes in weight post.df3$weight.change <- post.df3$weight.after - post.df3$weight post.swapped <- post.df3[post.df3$swapped == 1,] post.not.swapped <- post.df3[post.df3$swapped == 0,] # add OLRE post.df3$olre <- seq(1:nrow(post.df3)) ### Repeatabilites of behaviour and web structure ### ## Repeatability of behaviour ## attack.rpt <- rptBinary(attack ~ (1|ID) + weight.c + day + obs.day + distance.n + as.factor(web), grname = "ID", data = behav.df2, nboot = 1000, npermut = 1000, adjusted = TRUE) attack.rpt ## Repeatability of number of gumfooted lines ## gum.rpt <- rptPoisson(gum ~ (1|ID) + weight.c + as.factor(web), grname = "ID", data = web.df, nboot = 1000, npermut = 1000, adjusted = TRUE) gum.rpt ## Repeatability of structural lines ## str.rpt <- rptPoisson(structural ~ (1|ID) + weight.c + as.factor(web), grname = "ID", data = web.df, nboot = 1000, npermut = 1000, adjusted = TRUE) str.rpt ### Phenotypic trait correlations ### ## Relationship between GF and behaviour ## gum.attack <- glmer(gumfooted.own ~ (1|ID) + (1|web) + (1|olre) + weight.c + attack.number, data = meso.all.1, family = "poisson", control=glmerControl(optCtrl=list(maxfun = 100000), optimizer="bobyqa")) gum.attack.ci <- confint(gum.attack, method = "boot") # calculate CI's summary(gum.attack) # model summary gum.attack.ci[6:6,] # CI's for attack.number # Figure 3A par(mfrow = c(1,1)) plot(jitter(meso.all.1$gumfooted.own) ~ jitter(meso.all.1$attack.number), xlab = "Attack number", ylab = "Number of gumfooted lines", xlim = c(-0.5, 9.5), ylim = c(-1,85), cex.lab = 1.5, cex.axis = 1.5, cex = 1.5 ) seq.attack <- seq(min(meso.all.1$attack.number), max(meso.all.1$attack.number), 0.1) gf.att.coef <- as.data.frame(coef(summary(gum.attack))) # mean attack points(seq.attack, exp(gf.att.coef[1,1] + gf.att.coef[2,1]*0 + gf.att.coef[3,1]*seq.attack), type = "l", col = "black", lwd = 3) ## Relationship between SL and behaviour ## str.attack <- glmer(structural.own ~ (1|ID) + (1|web) + (1|olre) + weight.c + attack.number, data = meso.all.1, family = "poisson", control=glmerControl(optCtrl=list(maxfun = 100000), optimizer="bobyqa")) str.attack.ci <- confint(str.attack, method = "boot") # calculate CI's summary(str.attack) # model summary # plot it. This figure is not in the manuscript. plot(jitter(meso.all.1$structural.own) ~ jitter(meso.all.1$attack.number), xlab = "Attack number", ylab="Number of structural lines", xlim = c(-0.5, 9.5), cex.lab = 1.5, cex.axis = 1.5, cex = 1.5) seq.attack <- seq(min(meso.all.1$attack.number), max(meso.all.1$attack.number), 0.1) str.att.coef <- as.data.frame(coef(summary(str.attack))) # mean attack points(seq.attack, exp(str.att.coef[1,1] + str.att.coef[2,1]*0 + str.att.coef[3,1]*seq.attack) , type = "l", col = "black", lwd = 3) ## Relationship between GF and SL ## gum.str <- glmer(gumfooted.own ~ (1|ID) + (1|web) + (1|olre) + weight.c + structural.own, data = meso.all.1, family = "poisson", control=glmerControl(optCtrl=list(maxfun = 100000), optimizer="bobyqa")) gum.str.ci <- confint(gum.str, method = "boot") # calculate CI's summary(gum.str) # model summary # Figure 3B plot(jitter(meso.all.1$gumfooted.own) ~ jitter(meso.all.1$structural.own), xlab = "Number of structural lines", ylab="Number of gumfooted lines", ylim = c(-1,85), cex.lab = 1.5, cex.axis = 1.5, cex = 1.5) seq.str <- seq(min(meso.all.1$structural.own), max(meso.all.1$structural.own), 0.1) gf.str.coef <- as.data.frame(coef(summary(gum.str))) points(seq.str, exp(gf.str.coef[1,1] + gf.str.coef[2,1]*0 + gf.str.coef[3,1]*seq.str) , type = "l", col = "black", lwd = 3) ### Effect of transplantation on changes in web structure and body mass ### ## Change in gumfooted lines ## gum.change <- lmer(gumfooted.change ~ (1|ID) + weight.c + swapped, data = post.df3) # this model produces singular fit coefs.gum.change <- data.frame(coef(summary(gum.change))) coefs.gum.change$p.z <- 2 * (1 - pnorm(abs(coefs.gum.change$t.value))) coefs.gum.change # coefficients table gum.change.ci <- confint(gum.change, method = "boot") gum.change.ci[5:5, ] # CI's for effect of being swapped r.squaredGLMM((gum.change)) ## Change in structural lines ## str.change <- lmer(structural.change ~ (1|ID) + weight.c + swapped, data = post.df3) coefs.str.change <- data.frame(coef(summary(str.change))) coefs.str.change$p.z <- 2 * (1 - pnorm(abs(coefs.str.change$t.value))) coefs.str.change # coefficients table str.change.ci <- confint(str.change, method = "boot") str.change.ci[5:5, ] # CI's for effect of being swapped r.squaredGLMM(str.change) ## Change in weight ## weight.change <- lmer(weight.change ~ (1|ID) + weight.c + swapped, data = post.df3) coefs.weight.change <- data.frame(coef(summary(weight.change))) coefs.weight.change$p.z <- 2 * (1 - pnorm(abs(coefs.weight.change$t.value))) coefs.weight.change # coefficients table weight.change.ci <- confint(weight.change, method = "boot") weight.change.ci[5:5, ] # CI's for effect of being swapped r.squaredGLMM(weight.change) ### Variables influencing foraging success ### # specify model global.model3 <- glmer(prey.por ~ (1|ID) + (1|web.ID) + (1|web) + (1|olre) + weight.c + swapped + time.c + attack.number.c + gumfooted.c + I(attack.number.c^2) + I(gumfooted.c^2) + attack.number.c*gumfooted.c + swapped*gumfooted.c + swapped*attack.number.c + attack.number.c*gumfooted.c*swapped, data = meso.all, family = "binomial", weights = prey.number, na.action = na.pass, control = glmerControl(optCtrl=list(maxfun = 100000), optimizer = "bobyqa") ) summary(global.model3) # model summary global.model3.ci <- confint(global.model3, method = "boot") global.model3.ci # CI's r.squaredGLMM(global.model3) # run uncentered model for Figure S1A # rerun model, this time with uncentered GF global.model3.nc1 <- glmer(prey.por ~ (1|ID) + (1|web.ID) + (1|web) + (1|olre) + weight.c + swapped + time.c + attack.number.c + gumfooted + I(attack.number.c^2) + I(gumfooted^2) + attack.number.c*gumfooted + swapped*gumfooted + swapped*attack.number.c + attack.number.c*gumfooted*swapped, data = meso.all, family = "binomial", weights = prey.number, na.action = na.pass, control = glmerControl(optCtrl=list(maxfun = 100000), optimizer = "bobyqa") ) # run uncentered model for Figure - Figure S1B # rerun model, this time with uncentered attack.number global.model3.nc2 <- glmer(prey.por ~ (1|ID) + (1|web.ID) + (1|web) + (1|olre) + weight.c + swapped + time.c + attack.number + gumfooted.c + I(attack.number^2) + I(gumfooted.c^2) + attack.number*gumfooted.c + swapped*gumfooted.c + swapped*attack.number + attack.number*gumfooted.c*swapped, data = meso.all, family = "binomial", weights = prey.number, na.action = na.pass, control = glmerControl(optCtrl=list(maxfun = 100000), optimizer = "bobyqa") ) # Figure S1A plot(jitter(meso.all$prey.por) ~ jitter(meso.all$gumfooted), xlab = "Number of gumfooted lines", ylab="Proportion of ants eaten in interval", ylim = c(-0.05,1.05), cex.lab = 1.5, cex.axis = 1.5, cex = 1.5, col = rgb(0, 0, 0, alpha = 0.5) ) seq.gf <- seq(min(meso.all$gumfooted), max(meso.all$gumfooted), 0.1) capture.all.gum.coef1 <- as.data.frame(coef(summary(global.model3.nc1))) points(seq.gf, plogis(capture.all.gum.coef1[1,1] + # intercept capture.all.gum.coef1[2,1]*0 + # weight.c capture.all.gum.coef1[3,1]*0.5 + # swapped capture.all.gum.coef1[4,1]*0 + # time.c capture.all.gum.coef1[5,1]*0 + # attack.number.c capture.all.gum.coef1[6,1]*seq.gf + # GF capture.all.gum.coef1[7,1]*0 + # attack.number.c^2 capture.all.gum.coef1[8,1]*seq.gf + # GF^2 capture.all.gum.coef1[9,1]*0*seq.gf + # attack.number.c*GF capture.all.gum.coef1[10,1]*0.5*seq.gf + # swapped*GF capture.all.gum.coef1[11,1]*0.5*0 + # swapped*attack.number.c capture.all.gum.coef1[12,1]*0.5*0*seq.gf # three-way interaction ), type = "l", col = "black", lwd = 3) # Figure S1B plot(jitter(meso.all$prey.por) ~ jitter(meso.all$attack.number), xlab = "Foraging aggresssion", ylab="Proportion of ants eaten in interval", ylim = c(-0.05, 1.05), cex.lab = 1.5, cex.axis = 1.5, cex = 1.5, col = rgb(0, 0, 0, alpha = 0.5)) seq.attack <- seq(min(meso.all$attack.number), max(meso.all$attack.number), 0.1) capture.all.gum.coef2 <- as.data.frame(coef(summary(global.model3.nc2))) points(seq.attack, plogis(capture.all.gum.coef2[1,1] + # intercept capture.all.gum.coef2[2,1]*0 + # weight.c capture.all.gum.coef2[3,1]*0.5 + # swapped capture.all.gum.coef2[4,1]*0 + # time.c capture.all.gum.coef2[5,1]*seq.attack + # attack.number capture.all.gum.coef2[6,1]*0 + # GF.c capture.all.gum.coef2[7,1]*seq.attack + # attack.number^2 capture.all.gum.coef2[8,1]*0 + # GF.c^2 capture.all.gum.coef2[9,1]*seq.attack*0 + # attack.number*GF.c capture.all.gum.coef2[10,1]*0.5*0 + # swapped*GF.c capture.all.gum.coef2[11,1]*0.5*seq.attack + # swapped*attack.number capture.all.gum.coef2[12,1]*0.5*seq.attack*0 # three-way interaction ), type = "l", col = "black", lwd = 3) # Figure S2 means <- tapply(meso.all$interval.ants.eaten, as.factor(meso.all$swapped), mean) # calculate means medians <- tapply(meso.all$interval.ants.eaten, as.factor(meso.all$swapped), median) # calculate medians meso.all$swapped.named <- ifelse(meso.all$swapped == 1, 'Yes', 'No') stripchart(jitter(meso.all$interval.ants.eaten) ~ meso.all$swapped.named, vertical = TRUE, pch=1, cex = 0.5, cex.lab = 1.5, method="jitter", xlab = "Swapped", ylab = "Ants eaten in interval") points(x = 1:2,y = means, pch = 18, cex = 3, col = "blue") # add means points(x = 1:2,y = medians,pch = 19, cex = 3, col = "red") # add medians # run uncentered model for Figure S3 # rerun model, this time with uncentered time global.model3.nc3 <- glmer(prey.por ~ (1|ID) + (1|web.ID) + (1|web) + (1|olre) + weight.c + swapped + time + attack.number.c + gumfooted.c + I(attack.number.c^2) + I(gumfooted.c^2) + attack.number.c*gumfooted.c + swapped*gumfooted.c + swapped*attack.number.c + attack.number.c*gumfooted.c*swapped, data = meso.all, family = "binomial", weights = prey.number, na.action = na.pass, control = glmerControl(optCtrl=list(maxfun = 100000), optimizer = "bobyqa") ) # Figure S3 plot(jitter(meso.all$prey.por) ~ jitter(meso.all$time), xlab = "Time (hours)", ylab="Proportion of ants eaten in interval", ylim = c(-0.05, 1.05), cex.lab = 1.5, cex.axis = 1.5, cex = 1.5, col = rgb(0, 0, 0, alpha = 0.5) ) seq.time <- seq(min(meso.all$time), max(meso.all$time), 0.1) capture.all.gum.coef3 <- as.data.frame(coef(summary(global.model3.nc3))) points(seq.time, plogis(capture.all.gum.coef3[1,1] + # intercept capture.all.gum.coef3[2,1]*0 + # weight.c capture.all.gum.coef3[3,1]*0.5 + # swapped capture.all.gum.coef3[4,1]*seq.time + # time capture.all.gum.coef3[5,1]*0 + # attack.number.c capture.all.gum.coef3[6,1]*0 + # GF.c capture.all.gum.coef3[7,1]*0 + # attack.number.c*^2 capture.all.gum.coef3[8,1]*0 + # GF.c^2 capture.all.gum.coef3[9,1]*0*0 + # attack.number.c*GF.c capture.all.gum.coef3[10,1]*0.5*0 + # swapped*GF.c capture.all.gum.coef3[11,1]*0.5*0 + # swapped*attack.number.c capture.all.gum.coef3[12,1]*0.5*0*0 # three-way interaction ), type = "l", col = "black", lwd = 3) # Figure 4 pp <- cbind(meso.all$attack.number.c, meso.all$gumfooted.c, meso.all$prey.por) x1 <- pp[,1] # attack num x2 <- pp[,2] # GF y <- pp[,3] # prey.por x1seq = seq(-2.5, 1, 0.1) # attack.number sequence x2seq = seq(-1, 3.5, 0.1) # GF sequence f = function(x, y) {fixef(global.model3)[1] + # intercept x*fixef(global.model3)[5] + # attack.number.c y*fixef(global.model3)[6] + # GF.c x*fixef(global.model3)[7] + # attack.number.c^2 y*fixef(global.model3)[8] + # GF.c^2 x*y*fixef(global.model3)[9] # attack.number.c*GF.c } z = outer(x1seq, x2seq, f) z = exp(z)/(1 + exp(z)) res = persp(x1seq, x2seq, z, xlab = 'Foraging aggression (centered)', ylab = 'Gumfooted lines (centered)', zlab = 'Proportion of ants eaten', cex.lab = 1, theta = 40, phi = 15, r = 50, axes = T, cex.axis = 0.8, ticktype = 'detailed', zlim = c(-0.1, 1) ) # add points mypoints <- trans3d(x1, x2, y, pmat = res) points(mypoints, pch = 16, col = rgb(0,0,0,alpha = 0.5)) # add lines vpred <- trans3d(x1, x2, exp(f(x1,x2))/(1 + exp(f(x1,x2))), pmat = res) vobs <- trans3d(x1, x2, y, pmat = res) segments(vpred$x, vpred$y, vobs$x, vobs$y, col = "blue", lty = 3) ## AIC ## all.models <- dredge(global.model3) # this takes a while to run subset.models <- subset(all.models, delta <=2 ) # show only models with deltaAIC less than 2 summary(model.avg(subset.models)) # average only top models importance(subset.models)