##Analyses of D. ebraccatus tadpole phenotypic plasticity ##This script runs custom functions to resample tadpoles in order to calculate the magnitude of plasticity in response ##Instead of the bivariate model being run on different aspects of the phenotype, it is run within each measure but for ##responses to fish vs insect predators #Load packages library(lme4) library(MCMCglmm) library(car) library(tidyr) library(dplyr) library(broom) library(matrixStats) #Read in data TadPlast<-read.csv("TadPlast.csv") TadPlast<-na.omit(TadPlast) str(TadPlast) TadPlast$Treatment<-droplevels(TadPlast$Treatment) TadPlast$Treatment<-relevel(TadPlast$Treatment, ref="Fish") TadPlast$Family<-as.factor(TadPlast$Family) TadPlast$Pond<-droplevels(TadPlast$Pond) TadPlast<-droplevels(TadPlast) ################################################################################################################################################ ##Define functions to sample random sets of tadpoles and calculate the magnitude of plasticity (difference from the control) D.samp<-function(family, var) { temp.data<-TadPlast[TadPlast$Family==family,] temp.col<-which(names(temp.data)==var) D.temp<-temp.data[temp.data$TadNumber==sample(temp.data$TadNumber[temp.data$Treatment=="Insect"],1) & temp.data$Treatment=="Insect",temp.col] C.temp<-temp.data[temp.data$TadNumber==sample(temp.data$TadNumber[temp.data$Treatment=="Control"],1) & temp.data$Treatment=="Control",temp.col] D.plast<-D.temp-C.temp return(D.plast) } F.samp<-function(family, var) { temp.data<-TadPlast[TadPlast$Family==family,] temp.col<-which(names(temp.data)==var) F.temp<-temp.data[temp.data$TadNumber==sample(temp.data$TadNumber[temp.data$Treatment=="Fish"],1) & temp.data$Treatment=="Fish",temp.col] C.temp<-temp.data[temp.data$TadNumber==sample(temp.data$TadNumber[temp.data$Treatment=="Control"],1) & temp.data$Treatment=="Control",temp.col] F.plast<-C.temp-F.temp return(F.plast) } ##Resample the bivariate model for different phenotypic traits in Fish and Insect environments bi.mod.resamp<-function(families, var, reps) { output<-data.frame(matrix(data=NA, nrow=reps, ncol=8)) names(output)<-c("Run","Var","D.plast.herit","F.plast.herit","Mean.Gen.Corr","Mode.Gen.Corr","Gen.Corr.upper","Gen.Corr.lower") for(i in 1:reps) { #Build the dataset temp<-data.frame(matrix(nrow=340, ncol=6)) names(temp)<-c("animal","Family","FATHER","MOTHER","D.plast","F.plast") temp$animal<-1:340 temp$Family<-rep(unique(TadPlast$Family), each=10) temp$FATHER<-rep(unique(TadPlast$FATHER), each=10) temp$MOTHER<-rep(unique(TadPlast$MOTHER), each=20) for(j in 1:nrow(temp)) { temp$D.plast[j]<-D.samp(temp$Family[j], var) temp$F.plast[j]<-F.samp(temp$Family[j], var) } #Build the pedigree temp.ped<-data.frame(matrix(data=NA, nrow=391, ncol=3)) names(temp.ped)<-c("ID","MOTHER","FATHER") temp.ped[52:391,1]<-as.character(1:340) temp.ped[52:391,2]<-as.character(temp$MOTHER) temp.ped[52:391,3]<-as.character(temp$FATHER) temp.ped[1:51,1]<-as.character(TadPlast.ped[1:51,1]) #Run the model p.var<-matrix(c(var(temp$D.plast,na.rm=TRUE),0,0,var(temp$F.plast,na.rm=TRUE)),2,2) prior.bi1<-list(G=list(G1=list(V=p.var/2,n=2)),R=list(V=p.var/2,n=2)) model.bi1<-MCMCglmm(cbind(D.plast, F.plast) ~ trait-1, random = ~us(trait):animal, rcov=~us(trait):units, family = c("gaussian", "gaussian"), pedigree = temp.ped, data = temp, prior = prior.bi1, nitt=250000, thin=100, burnin= 50000, verbose = F) #Calculate heritability of each trait and the genetic correlation herit1<-model.bi1$VCV[,'traitD.plast:traitD.plast.animal']/(model.bi1$VCV[,'traitD.plast:traitD.plast.animal']+ model.bi1$VCV[,'traitD.plast:traitD.plast.units']) herit2<-model.bi1$VCV[,'traitF.plast:traitF.plast.animal']/(model.bi1$VCV[,'traitF.plast:traitF.plast.animal']+ model.bi1$VCV[,'traitF.plast:traitF.plast.units']) corr.gen<-model.bi1$VCV[,'traitD.plast:traitF.plast.animal']/sqrt(model.bi1$VCV[,'traitD.plast:traitD.plast.animal']* model.bi1$VCV[,'traitF.plast:traitF.plast.animal']) #Put it in the output output[i,1]<-i output[i,2]<-as.character(var) output[i,3]<-mean(herit1) output[i,4]<-mean(herit2) output[i,5]<-mean(corr.gen) output[i,6]<-posterior.mode(corr.gen) output[i,7]<-HPDinterval(corr.gen)[1] output[i,8]<-HPDinterval(corr.gen)[2] } return(output) } Saturation.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="Saturation", reps=500) write.csv(Saturation.bi.resamp, "Saturation.bi.resamp.csv") MaxTailDepth.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="MaxTailDepth", reps=500) write.csv(MaxTailDepth.bi.resamp, "MaxTailDepth.bi.resamp.csv") TailMuscleHeight.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="TailMuscleHeight", reps=500) write.csv(TailMuscleHeight.bi.resamp, "TailMuscleHeight.bi.resamp.csv") Hue.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="Hue", reps=500) write.csv(Hue.bi.resamp, "Hue.bi.resamp") BodyLength.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="BodyLength", reps=500) write.csv(BodyLength.bi.resamp, "BodyLength.bi.resamp.csv") TailMuscleWidth.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="TailMuscleWidth", reps=500) write.csv(TailMuscleWidth.bi.resamp, "TailMuscleWidth.bi.resamp.csv") TailLength.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="TailLength", reps=500) write.csv(TailLength.bi.resamp, "TailLength.bi.resamp.csv") TailSpotArea.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="TailSpotArea", reps=500) write.csv(TailSpotArea.bi.resamp, "TailSpotArea.bi.resamp.csv") HeadWidth.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="HeadWidth", reps=500) write.csv(HeadWidth.bi.resamp, "HeadWidth.bi.resamp.csv") TotalLength.bi.resamp<-bi.mod.resamp(families=unique(TadPlast$Family), var="TotalLength", reps=500) write.csv(TotalLength.bi.resamp, "TotalLength.bi.resamp.csv")