### Code for Paniw et al. "Demographic traits improve predictions of spatiotemporal changes in community resilience to drought" # R version 4.0.0 # date: Nov 15, 2020 # This code plots the three resilience mearures as a function of individual traits rm(list=ls(all=TRUE)) library(ggplot2) library(ggrepel) library(dplyr) library(stringr) library(Epi) library(vegan) library(goeveg) library(viridis) # Species cover per plot and year data=read.csv("/Users/maria/Dropbox/demoDonana/SuppMat/species_cover.csv") # Total cover per plot and year and initial impact (dieOffCont) cat=read.csv("/Users/maria/Dropbox/demoDonana/SuppMat/plot_cover.csv") colnames(cat)[1]="site" colnames(cat)[3]="dieOffCont" data$site=sapply(str_split(data$Plot,"_"), "[", 1) data$cond=sapply(str_split(data$Plot,"_"), "[", 2) data$dieoffCat=left_join(data,cat,by="site")$dieoffCat data$dieoffCont=left_join(data,cat,by="site")$dieOffCont #### DISSIMILARITY (Bray Curtis) #transform data mydata=decostand(data[,2:12],"max") row.names(mydata)=data$Plot diss=vegdist(mydata, method="bray", binary=FALSE, diag=F, upper=T, na.rm = FALSE) head(diss) # ORDINATION: Nonmetric multidimensional scaling # the lower the stress the better out=dimcheckMDS(mydata, k = 9, try = 100, autotransform = TRUE)# three is a good number of K comp=metaMDS(mydata,k=3,try = 20, trymax = 100) stressplot(comp) plot(comp, type='n') points(comp, display=c('sites'), choices=c(1, 2), pch=3, col="red") text(comp, display=c('species'), choices=c(1, 2), col='blue', cex=0.7) sampleScores <- comp$points # NMDS space distance dist2008=diag(as.matrix(dist(sampleScores,diag=T))[37:54,1:18]) dist2010=diag(as.matrix(dist(sampleScores,diag=T))[55:72,1:18]) dist2013=diag(as.matrix(dist(sampleScores,diag=T))[73:90,1:18]) dist2017=diag(as.matrix(dist(sampleScores,diag=T))[91:108,1:18]) dist2019=diag(as.matrix(dist(sampleScores,diag=T))[109:dim(as.matrix(diss))[1],1:18]) # TRAIT DATA demo=read.csv("/Users/maria/Dropbox/demoDonana/SuppMat/demoTraits.csv") demo=demo[c(5,3,1,4,2,10,11,6,9,8,7),] # add minimum size at flowering df=read.csv("/Users/maria/Dropbox/demoDonana/SuppMat/seed_flower.csv") head(df) df$Flower=factor(df$Flower) levels(df$Flower)=c("0","1") df$Flower=as.numeric(as.character(df$Flower)) df$size=log(df$Diameter.cm.*df$h..cm.) hist(df$size) # Loop through species min.size.fl=data.frame(Species=rep(NA,length(unique(df$species))), min.size.obs=rep(NA,length(unique(df$species))), min.size.ROC=rep(NA,length(unique(df$species)))) for(i in 1: length(unique(df$species))){ df2=df[df$species==unique(df$species)[i],] min.size.fl$Species[i]=as.character(unique(df$species)[i]) rc <- ROC(form = df2$Flower ~ df2$size, plot="sp") ## optimal combination opt <- which.max(rowSums(rc$res[, c("sens", "spec")])) # what size min.size.fl$min.size.ROC[i]= df2$size[which(predict(rc$lr,type="response")==rc$res$lr.eta[opt])][1] # observed mimimum min.size.fl$min.size.obs[i]=min(df2$size[df2$Flower==1]) } min.size.fl=min.size.fl[order(min.size.fl$min.size.ROC),] min.size.fl$Species=factor(min.size.fl$Species,levels=min.size.fl$Species[order(min.size.fl$min.size.ROC)]) levels(min.size.fl$Species)=c("Helycrissum picardii", "Halimium commutatum","Lavandula stoechas","Thymus mastichina", "Cistus libanotis","Erica scoparia","Staurachantus genistoides","Rosmarinus officinalis", "Halimium halimifolium","Ulex australis","Juniperus phoenicea") demo$minSizeAtFlower=left_join(demo,min.size.fl,by="Species")$min.size.ROC cover=data[,2:12] ### Community weighted means ageFR=NULL long1=NULL sizeFL=NULL r_a_ratio=NULL for(i in 1:nrow(cover)){ ageFR[i]=sum(cover[i,]*demo[,2])/sum(cover[i,]) long1[i]=sum(cover[i,]*demo[,3])/sum(cover[i,]) sizeFL[i]=sum(cover[i,]*demo[,5])/sum(cover[i,]) r_a_ratio[i]=sum(cover[i,]*demo[,4])/sum(cover[i,]) } traits=cbind(ageFR, long1, sizeFL, r_a_ratio) # ADD FUNCTIONAL TRAITS ft=read.csv("/Users/maria/Dropbox/demoDonana/SuppMat/funcTraits.csv") ft=ft[c(5,3,1,4,2,10,11,6,9,8,7),] LDMC=NULL leaf_chlor=NULL SLA=NULL N_leaf=NULL C13=NULL SRA_root=NULL RDMC=NULL seed.size=NULL height=NULL dry_matter=NULL for(i in 1:nrow(cover)){ LDMC[i]=sum(cover[i,]*ft[,2])/sum(cover[i,]) leaf_chlor[i]=sum(cover[i,]*ft[,3])/sum(cover[i,]) SLA[i]=sum(cover[i,]*ft[,4])/sum(cover[i,]) N_leaf[i]=sum(cover[i,]*ft[,5])/sum(cover[i,]) C13[i]=sum(cover[i,]*ft[,6])/sum(cover[i,]) SRA_root[i]=sum(cover[i,]*ft[,7])/sum(cover[i,]) RDMC[i]=sum(cover[i,]*ft[,8])/sum(cover[i,]) seed.size[i]=sum(cover[i,]*ft[,9])/sum(cover[i,]) height[i]=sum(cover[i,]*ft[,10])/sum(cover[i,]) dry_matter[i]=sum(cover[i,]*ft[,11])/sum(cover[i,]) } traits=cbind(traits,LDMC,leaf_chlor,SLA,N_leaf, C13,SRA_root,RDMC,seed.size,height,dry_matter) traits[is.na(traits)]=0 options(na.action = "na.fail") ## RELATIVE RESILIENCE IN TOTAL COVER res2008=(cat$cover.2008-cat$cover.at.2007)/cat$cover.pre.2007 res2010=(cat$cover.2010-cat$cover.at.2007)/cat$cover.pre.2007 res2013=(cat$cover.2013-cat$cover.at.2007)/cat$cover.pre.2007 res2017=(cat$cover.2017-cat$cover.at.2007)/cat$cover.pre.2007 res2019=(cat$cover.2019-cat$cover.at.2007)/cat$cover.pre.2007 # put data together for models new.dat=data.frame(RRs=c(res2008,res2010,res2013,res2017,res2019), time=rep(c(1,3,6,10,12),each=18), site=rep(substr(data$site[1:18],1,2),5), siteAff=rep(data$site[1:18],5)) ## ABSOLUTE RESILIENCE IN TOTAL COVER res2008=cat$cover.2008/cat$cover.pre.2007 res2010=cat$cover.2010/cat$cover.pre.2007 res2013=cat$cover.2013/cat$cover.pre.2007 res2017=cat$cover.2017/cat$cover.pre.2007 res2019=cat$cover.2019/cat$cover.pre.2007 new.dat$Rs=c(res2008,res2010,res2013,res2017,res2019) ## RESILIENCE IN COMPOSITION Rt=diag(as.matrix(dist(sampleScores,diag=T))[19:36,1:18])# distance pre.drought - drought res2008=dist2008 res2010=dist2010 res2013=dist2013 res2017=dist2017 res2019=dist2019 new.dat$RRc=c(1-(res2008/Rt),1-(res2010/Rt),1-(res2013/Rt),1-(res2017/Rt),1-(res2019/Rt)) # Add traits new.dat=cbind(new.dat,traits[c(37:126),]) new.dat$site=factor(new.dat$site) ### Plots # Put the data into long format the brtue way df=new.dat[,1:7] df$trait=colnames(df)[7] colnames(df)[7]="value" for(i in 8:ncol(new.dat)){ temp=new.dat[,c(1:6,i)] temp$trait=colnames(new.dat)[i] colnames(temp)[7]="value" df=rbind(df,temp) } df$aff[df$siteAff%in%c("CM1","CM2","OM1","OM4","MM1","MM5")]="severe" df$aff[df$siteAff%in%c("CM4","CM5","OM3","OM6","MM3","MM6")]="medium" df$aff[df$siteAff%in%c("CM3","CM6","OM2","OM5","MM2","MM4")]="low" df$Year=factor(df$time) levels(df$Year)=c(2008,2010,2013,2017,2019) df$trait=factor(df$trait,levels=colnames(new.dat)[c(7:18,20,19)]) levels(df$trait)=c("Age[FR]","L","Size[FR]","R/A","LDMC","LChl","SLA","LNC",expression(paste(delta^{13},"C")),"SRA","RDMC","S[mass]","SDMC","Phg") # Absolute resilience a=ggplot(df, aes(value,Rs,col=aff))+ geom_point(data=df, aes(value,Rs,shape=Year),size=2)+ facet_wrap(.~trait,ncol=4,scales = "free_x",labeller = label_parsed)+ geom_smooth(data=df, aes(value,Rs),method = "lm", formula = y ~ x, se = T,inherit.aes = F)+ scale_color_viridis("Initial impact",discrete=T,option="C",begin=0,end=0.9)+ xlab("Trait value")+ ylab("Rs")+ theme_bw()+ theme(panel.grid = element_blank())+ theme(axis.text = element_text(size=20,colour="black"))+ theme(axis.title = element_text(size=24,colour="black"))+ theme(axis.title.x = element_text(vjust=-0.4), axis.title.y = element_text(vjust=2.3))+ theme(legend.title = element_text(size=20), legend.text = element_text(size=18), legend.position = "top")+ theme(plot.margin = unit(c(1,0.5,0.5,0.5), "cm"))+ theme(strip.text = element_text(size = 20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background = element_blank(), plot.title = element_text(size=22)) a # Relative resilience a=ggplot(df, aes(value,RRs,col=aff))+ geom_point(data=df, aes(value,RRs,shape=Year),size=2)+ facet_wrap(.~trait,ncol=4,scales = "free_x",labeller = label_parsed)+ geom_smooth(data=df, aes(value,RRs),method = "lm", formula = y ~ x, se = T,inherit.aes = F)+ scale_color_viridis("Initial impact",discrete=T,option="C",begin=0,end=0.9)+ xlab("Trait value")+ ylab("RRs")+ theme_bw()+ theme(panel.grid = element_blank())+ theme(axis.text = element_text(size=20,colour="black"))+ theme(axis.title = element_text(size=24,colour="black"))+ theme(axis.title.x = element_text(vjust=-0.4), axis.title.y = element_text(vjust=2.3))+ theme(legend.title = element_text(size=20), legend.text = element_text(size=18), legend.position = "top")+ theme(plot.margin = unit(c(1,0.5,0.5,0.5), "cm"))+ theme(strip.text = element_text(size = 20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background = element_blank(), plot.title = element_text(size=22)) a # Resilience in composition a=ggplot(df, aes(value,RRc,col=aff))+ geom_point(data=df, aes(value,RRc,shape=Year),size=2)+ facet_wrap(.~trait,ncol=4,scales = "free_x",labeller = label_parsed)+ geom_smooth(data=df, aes(value,RRc),method = "lm", formula = y ~ x, se = T,inherit.aes = F)+ scale_color_viridis("Initial impact",discrete=T,option="C",begin=0,end=0.9)+ xlab("Trait value")+ ylab("RRc")+ theme_bw()+ theme(panel.grid = element_blank())+ theme(axis.text = element_text(size=20,colour="black"))+ theme(axis.title = element_text(size=24,colour="black"))+ theme(axis.title.x = element_text(vjust=-0.4), axis.title.y = element_text(vjust=2.3))+ theme(legend.title = element_text(size=20), legend.text = element_text(size=18), legend.position = "top")+ theme(plot.margin = unit(c(1,0.5,0.5,0.5), "cm"))+ theme(strip.text = element_text(size = 20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.background = element_blank(), plot.title = element_text(size=22)) a