######################################## # Code to reproduce the figures from "The ghost of disturbance past: long-term effects of pulse disturbances on community biomass and composition" # Authors: Jacquet, C., Altermatt, F. ############################################## ############################################## # Clear variable space rm(list=ls()) # Set working directory: #setwd("~/Downloads/doi_10.5061_dryad.zkh189378") # Load libraries and functions library(RColorBrewer) library(scico) library(ggplot2) library(grid) library(scales) library(sciplot) # For the heatmaps of figures 2-4: source("multiplot.R") source("param_heatmaps_biomass.R") source("param_heatmaps_sp.R") source("stats_fun.R") # Load data load("community_biomass.Rdata") # Data Frame with community biomass at each time point (day=1 to day=21, and day=39) for all microcosms. PA = read.table("table_species_PA.txt",sep="\t",header=T) # Data Frame with species presence/absence at day=39 for all microcosms. p_value=0.1 # set p_value dat1=community_biomass[-which(community_biomass$day==39),] # dat1 is community biomass between day=1 and day=21. dat2=community_biomass[which(community_biomass$day==39),] # dat2 is community biomass at day=39. # Graph specifications int = levels(community_biomass$int)[1:5] freq = levels(community_biomass$freq)[2:5] treat=type.convert(expand.grid(int,freq)) colnames(treat)=c("int","freq") number=as.factor(round(c(1/3,1/6,1/9,1/12),2)) coul=c("RdBu") # color gradient titles = c("Direct effect","Legacy effect") tags =c("(a)","(b)") my.rect=cbind(expand.grid(seq(0.5,4.5,1),seq(0.5,3.5,1)),expand.grid(seq(0.5,4.5,1)+1,seq(0.5,3.5,1)+1)) colnames(my.rect)=c("xmin","ymin","xmax","ymax") color = rep(c("#d73027","#f46d43","#fdae61","#74add1","#4575b4"),4) color=color[20:1] color2=rep(c("#d7302750","#f46d4350","#fdae6150","#74add150","#4575b450"),4) color2=color2[20:1] ######################################################### # FIGURE 1 ######################################################### D_day=list() # Disturbance days for each frequency for(i in 1 :5){ D_day[[i]]=c(1,4,7,10,13,16,19) } for(i in 6:10){ D_day[[i]]=c(1,7,13,19) } for(i in 11:15){ D_day[[i]]=c(1,10,19) } for(i in 16:20){ D_day[[i]]=c(1,13) } BOUND=NULL # Y-axis limits for each treatment X0 = as.factor(community_biomass$day[which(community_biomass$freq == "Co")]) Y0 = community_biomass$bioarea_per_volume[which(community_biomass$freq == "Co")] for (j in 1:length(freq)) { for (i in 1:length(int)) { X = as.factor(community_biomass$day[which(community_biomass$freq == freq[j] & community_biomass$int == int[i])]) Y = community_biomass$bioarea_per_volume[which(community_biomass$freq == freq[j] & community_biomass$int == int[i])] pdf(file=NULL) graph = lineplot.CI(X,Y,file=NULL) graph2 = lineplot.CI(X0,Y0,file=NULL) dev.off() BOUND = rbind(BOUND,range(c(graph2$CI,graph$CI))) } } #png("Figure1.png", width = 8, height = 4, units = 'in', res = 300) par(mfrow=c(4,5),mar=rep(1.5,4),oma=c(3,3,0.5,0.5)) X0 = as.factor(community_biomass$day[which(community_biomass$freq == "Co")]) Y0 = community_biomass$bioarea_per_volume[which(community_biomass$freq == "Co")] for (i in 1:20){ X = as.factor(community_biomass$day[which(community_biomass$freq == as.character(treat$freq[i]) & community_biomass$int == as.character(treat$int[i]))]) Y = community_biomass$bioarea_per_volume[which(community_biomass$freq == as.character(treat$freq[i]) & community_biomass$int == as.character(treat$int[i]))] test=lineplot.CI(X0,Y0,err.width=0.03,pch=16,col="#63636380",err.col="#63636380",lty=1,lwd=1.5,xlab=NA,ylab=NA,err.lty=1,x.cont=T,ylim=BOUND[i,],cex.axis=1.3) abline(v=D_day[[i]],lty=2,col="#63636390") par(new=TRUE) lineplot.CI(X,Y,err.width=0.03,pch=16,col=color[i],err.col=color2[i],lty=1,lwd=1.5,xlab=NA,ylab=NA,err.lty=1,x.cont=T,ylim=BOUND[i,],xaxt="n",yaxt="n",bty="n") } mtext(side=1,at=0.5,outer=T,line=1.5,"Time (day)",cex=1.2) mtext(side=2,at=0.5,outer=T,line=0.8, expression("Community biomass " (mu*m^2/mu*l)),cex=1.2) #dev.off() ######################################################### # ######################################################### # FIGURE 2 ######################################################### dat = dat1 va = (dat$bioarea_per_volume) title = titles[1] tag = tags[1] Va_name = c("Biomass") res_press_biomass= stats_fun(dat,int,freq,va) # statistics for Table S2 p1 = param_heatmaps_biomass(dat,int,freq,number,va,p_value,coul,title,Va_name,tag) # panel (a) of Figure 2 dat = dat2 va = dat$bioarea_per_volume title = titles[2] tag = tags[2] Va_name = c("Biomass") p2 = param_heatmaps_biomass(dat,int,freq,number,va,p_value,coul,title,Va_name,tag) # statistics for Table S2 res_lt_biomass= stats_fun(dat,int,freq,va) # panel (b) of Figure 2 #png("Figure2.png", width = 4, height = 4, units = 'in', res = 300) multiplot(p1,p2,cols=1) #dev.off() ######################################################### ######################################################### # Figure 3 ######################################################### # 1) Statistics for Table S2 SR=apply(PA[,3:12],1,sum) # Species richness dat=dat2 va = SR Va_name = c("Species richness") X = va[which(dat$freq == "Co")] res_SR = stats_fun(dat,int,freq,va) res_SR[,3]=paste(" ",as.character(res_SR[,3])) res_SR[,4]=paste(" ",as.character(res_SR[,4])) res_SR[,6]=paste(" ",as.character(res_SR[,6])) # Figure 3 dat=dat2 va = SR Va_name = c("Species richness") X = va[which(dat$freq == "Co")] RES = NULL for (i in 1:length(int)) { for (j in 1:length(freq)) { Y = va[which(dat$freq == freq[j] & dat$int == int[i])] test=t.test(X,Y) res=cbind(freq[j],int[i],test$p.value,test$estimate[2]) RES=rbind(RES,res) }} data=data.frame(expand.grid(number,int),as.numeric(RES[,4])) colnames(data)=c("Number","Intensity",Va_name) p_val= my.rect[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value),] diff=as.numeric(RES[,4])-mean(X) couls = rep(" ",20) couls[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value & diff>0)]="+" couls[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value & diff<0)]="-" data2=data p3=ggplot(data2) + geom_raster(aes(y=Number,x=Intensity,fill=data2[,3]))+ scale_fill_gradientn(name=Va_name,colours=(brewer.pal(9,coul)))+ coord_fixed()+ geom_rect(data=my.rect,aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=0,color=1,lwd=0.2)+ geom_text(data=data,aes(y=Number,x=Intensity,label = couls),size=5)+ labs(x = "Disturbance intensity", y = "Disturbance frequency") # png("Figure3.png", width = 4, height = 4, units = 'in', res = 300) p3 # dev.off() ######################################################### ######################################################### # FIGURE 4 ######################################################### dat=dat2 va = PA$Ble title = expression(italic("Blepharisma sp.")) #c("Blepharisma sp.") tag= "(a)" Va_name = c("P(X=1)") p1 = param_heatmaps_sp(dat,int,freq,number,va,p_value,coul,title,Va_name,tag,c(" "),c("Disturbance frequency")) #p1 ########################### va = PA$Cep..rotifer. title = expression(italic("Cephalodella sp.") ) tag= "(b)" Va_name = c("P(X=1)") p2 =param_heatmaps_sp(dat,int,freq,number,va,p_value,coul,title,Va_name,tag,c(" "),c(" ")) #p2 ########################### va = PA$Chi title = expression(italic("Chilomonas sp.") ) tag= "(c)" Va_name = c("P(X=1)") p3 =param_heatmaps_sp(dat,int,freq,number,va,p_value,coul,title,Va_name,tag,c(" "),c(" ")) #p3 ########################### va = PA$Col.Pau.Pca title = expression(italic("Colpidium + Paramecium") ) Va_name = c("P(X=1)") tag= "(d)" p4 = param_heatmaps_sp(dat,int,freq,number,va,p_value,coul,title,Va_name,tag,c("Disturbance intensity"),c("Disturbance frequency")) #p4 ########################### ########################### va = PA$Eug title = expression(italic("Euglena gracilis ") ) Va_name = c("P(X=1)") tag= "(e)" X = va[which(dat$freq == "Co")] RES = NULL for (i in 1:length(int)) { for (j in 1:length(freq)) { Y = va[which(dat$freq == freq[j] & dat$int == int[i])] if (mean(Y)!=1){ test=t.test(Y,mu=X[1]) res=cbind(freq[j],int[i],test$p.value,test$estimate)} else {res=cbind(freq[j],int[i],1,1)} RES=rbind(RES,res) } } data=data.frame(expand.grid(number,int),as.numeric(RES[,4])) colnames(data)=c("Number","Intensity",Va_name) p_val= my.rect[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value),] diff=as.numeric(RES[,4])-mean(X) couls = rep(" ",20) couls[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value & diff>0)]="+" couls[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value & diff<0)]="-" data_Eug=data p5=ggplot(data_Eug) + geom_raster(aes(y=Number,x=Intensity,fill=data_Eug[,3]))+ scale_fill_gradientn(name=" ",colours=(brewer.pal(9,coul)))+ coord_fixed()+ geom_rect(data=my.rect,aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=0,color=1,lwd=0.2)+ geom_text(data=data,aes(y=Number,x=Intensity,label = couls),size=5)+ labs(x = c("Disturbance intensity"), y = c("Disturbance frequency"),title=title,tag=tag) #p5 # ########################### va = PA$Eup title = expression(italic("Euplotes aediculatus")) Va_name = c("P(X=1)") tag= "(f)" X = va[which(dat$freq == "Co")] RES = NULL for (i in 1:length(int)) { for (j in 1:length(freq)) { Y = va[which(dat$freq == freq[j] & dat$int == int[i])] if (mean(Y)!=1){ test=t.test(Y,mu=X[1]) res=cbind(freq[j],int[i],test$p.value,test$estimate)} else {res=cbind(freq[j],int[i],1,1)} RES=rbind(RES,res) } } data=data.frame(expand.grid(number,int),as.numeric(RES[,4])) colnames(data)=c("Number","Intensity",Va_name) p_val= my.rect[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value),] diff=as.numeric(RES[,4])-mean(X) couls = rep(" ",20) couls[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value & diff>0)]="+" couls[which(p.adjust(as.numeric(RES[,3]),"fdr")<=p_value & diff<0)]="-" data_Eup=data p6=ggplot(data_Eup) + geom_raster(aes(y=Number,x=Intensity,fill=data_Eup[,3]))+ scale_fill_gradientn(name=" ",colours=(brewer.pal(9,coul)))+ coord_fixed()+ geom_rect(data=my.rect,aes(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax),alpha=0,color=1,lwd=0.2)+ geom_text(data=data,aes(y=Number,x=Intensity,label = couls),size=5)+ labs(x = c("Disturbance intensity"), y = c("Disturbance frequency"),title=title,tag=tag) #p6 # Figure 4 # png("Figure4.png", width = 4, height = 4, units = 'in', res = 300) multiplot(p1,p4,p2,p5,p3,p6,cols=3) # dev.off() ######################################################### ######################################################### # FIGURE 5 ######################################################### I=c(0.1,0.3,0.5,0.7,0.9) T=c(3,6,9,12) treatment=expand.grid(T,I) crit=-log(1-treatment[,2])/treatment[,1] coul=c("#67000d99","#cb181d99","#fb6a4a99","#fcbba199") coul2=c("#67000d","#cb181d","#fb6a4a","#fcbba1") int=seq(0,1,0.01) freq=c(3,6,9,12) treatment2=expand.grid(freq,int) crit2=-log(1-treatment2[,2])/treatment2[,1] # Figure 4 # png("Figure5.png", width = 4, height = 4, units = 'in', res = 300) plot(treatment$Var2*100,(crit),xlab="Disturbance intensity (%)",ylab="Disturbance regime strength",pch=NA,xlim=c(0,100),las=1,cex.lab=1.3,cex.axis=1.3) f=3 points(treatment$Var2[which(treatment$Var1==f)]*100,(crit[which(treatment$Var1==f)]),col=coul2[1],pch=19) lines(treatment2$Var2[which(treatment2$Var1==f)]*100,(crit2[which(treatment2$Var1==f)]),col=coul[1],lwd=1.5) f=6 points(treatment$Var2[which(treatment$Var1==f)]*100,(crit[which(treatment$Var1==f)]),col=coul2[2],pch=19) lines(treatment2$Var2[which(treatment2$Var1==f)]*100,(crit2[which(treatment2$Var1==f)]),col=coul[2],lwd=1.5) f=9 points(treatment$Var2[which(treatment$Var1==f)]*100,(crit[which(treatment$Var1==f)]),col=coul2[3],pch=19) lines(treatment2$Var2[which(treatment2$Var1==f)]*100,(crit2[which(treatment2$Var1==f)]),col=coul[3],lwd=1.5) f=12 points(treatment$Var2[which(treatment$Var1==f)]*100,(crit[which(treatment$Var1==f)]),col=coul2[4],pch=19) lines(treatment2$Var2[which(treatment2$Var1==f)]*100,(crit2[which(treatment2$Var1==f)]),col=coul[4],lwd=1.5) legend("topleft",c("0.33","0.17","0.11","0.08"),title="Disturbance frequency",col=coul2,bty="n",lty=1,lwd=1.3,cex=1.2) # dev.off() #########################################################