Note: This code is taken from an R Markdown Notebook. The preview shows you a rendered HTML copy of the contents of the editor.

Presets and reading data

Environmental data for the plots and an OTU table were read from disk (to run the notebook, contact AHB).

A presence/absence version of the OTU-table is prepared, as well as a rarefied version and a presence/absence version of this. In addition, OTU tables containing only OTUs of Glomus or of Claroidoglomus were prepared. Colours for the subplots are chosen from the “spectral” palette.

library(RColorBrewer)
library(vegan)
## Warning: package 'vegan' was built under R version 3.4.4
spcolsP <- brewer.pal(6,"Spectral")
spcols <- as.vector(t(unlist(sapply(spcolsP,function(x)colorRampPalette(c("grey90",x))(7)[3:7]))))
subplots <- paste("subplot",1:30,sep="_")
tps <- 1:6
tabPlots <- data.frame("sample"=1:360,
                       "tp"=rep(1:6,each=60),
                       "subplot"=paste("subplot",rep(rep(1:30,each=2),times=6),sep="_"),
                       stringsAsFactors=F)

env <- AMF[,1:58]
envSubplot <- rep(as.factor(rep(1:30,each=2)),6)
OTUS <- AMF[,59:213]

otuM <- as.matrix(otuPre[,4:ncol(otuPre)])
rownames(otuM) <- otuPre$SampleID
otuMR <- rrarefy(otuM,min(rowSums(otuM)))
otuB <- matrix(as.numeric(otuM>0),ncol=ncol(otuM),nrow=nrow(otuM),dimnames = dimnames(otuM))
otuRB <- matrix(as.numeric(otuMR>0),ncol=ncol(otuMR),nrow=nrow(otuMR),dimnames = dimnames(otuMR))

glomusM <- as.matrix(glomusPre[,4:ncol(glomusPre)])
rownames(glomusM) <- glomusPre$SampleID
glomusB <- matrix(as.numeric(glomusM>0),ncol=ncol(glomusM),nrow=nrow(glomusM),dimnames = dimnames(glomusM))

claroM <- as.matrix(claroPre[,4:ncol(claroPre)])
rownames(claroM) <- claroPre$SampleID
claroB <- matrix(as.numeric(claroM>0),ncol=ncol(claroM),nrow=nrow(claroM),dimnames = dimnames(claroM))

scheme <- list()
scheme$sites <- data.frame("x"=jitter(rep(1:6,times=5),amount=0.2),"y"=jitter(rep(5:1,each=6),amount=0.2),stringsAsFactors = F)

{plot(rep(1:6,times=5),rep(5:1,each=6),col=spcols,pch=16,cex=5,axes=F,ann=F,xlim=c(0.2,6.8),ylim=c(0.2,5.8))
text(rep(1:6,times=5),rep(5:1,each=6),labels=1:30)}

A few functions are custom-defined. In particular, we make use of the betapart package - however, the functions adjust the comparisons between samples to take into account the design which has always two samples taken at each time point in the same subplot.

#basic module:
subset.betapart.core <- function(bpc, rows) 
{
  bpc$data <- bpc$data[rows, ]
  bpc$shared <- bpc$shared[rows, rows]
  bpc$not.shared <- bpc$not.shared[rows, rows]
  bpc$sumSi <- sum(diag(bpc$shared))
  bpc$St <- sum(colSums(bpc$data) > 0)
  bpc$a <- bpc$sumSi - bpc$St
  bpc$sum.not.shared <- bpc$sum.not.shared[rows, rows]
  bpc$max.not.shared <- bpc$max.not.shared[rows, rows]
  bpc$min.not.shared <- bpc$min.not.shared[rows, rows]
  return(bpc)
}

#for p/a data
beta.sample.kezia <- function (x, index.family = "sorensen", sites = nrow(x$data), tps = 6, repl = 2, samples = 64) 
{
  if (!inherits(x, "betapart")) {
    x <- betapart.core(x)
  }
  if(samples > repl^tps) {
    #print(paste("Number of requested samples greater than maximum number of permutations. Max. permutations",
    #            repl^tps, "used."))
    samples <- repl^tps
  }
  #pb <- txtProgressBar(min = 0, max = samples, style = 3)
  if (sites > nrow(x$data)) 
    stop("More sites requested for sample than are in the dataset")
  index.family <- match.arg(index.family, c("jaccard", "sorensen"))
  results.n <- as.data.frame(matrix(nrow = samples, ncol = 3))
  if(repl == 2){
    sampleVec <- sample(1:(repl^tps),samples)
    for (i in 1:samples) {
      pos <- as.integer(intToBits(sampleVec[i]-1))[1:tps]
      sample.pos <- 2*1:tps-pos
      x.sample <- subset.betapart.core(x, sample.pos)
      x.beta <- beta.multi(x.sample, index.family)
      results.n[i, ] <- unlist(x.beta)
      #setTxtProgressBar(pb, i)
    }
  }else{
    stop("function is only implemented for two replicates")
  }  
  #close(pb)
  names(results.n) <- names(x.beta)
  result <- list(sampled.values = results.n, mean.values = sapply(results.n, 
                                                                  mean), sd.values = sapply(results.n, sd))
  return(result)
}

#twin for abundance data (not used in the notebook)
# beta.sample.abund.kezia <- function (x, index.family = "bray", tps = 6, repl = 2, samples = 64) 
# {
#   index.family <- match.arg(index.family, c("bray", "ruzicka"))
#   if(samples > repl^tps) {
#     # print(paste("Number of requested samples greater than maximum number of permutations. Max. permutations",
#     #              repl^tps, "used."))
#     samples <- repl^tps
#   }
#   #pb <- txtProgressBar(min = 0, max = samples, style = 3)
#   results.n <- as.data.frame(matrix(nrow = samples, ncol = 3))
#   if(repl == 2){
#     sampleVec <- sample(1:(repl^tps),samples)
#     for (i in 1:samples) {
#       pos <- as.integer(intToBits(sampleVec[i]-1))[1:tps]
#       sample.pos <- 2*1:tps-pos
#       x.sample <- x[sample.pos,]
#       x.beta <- beta.multi.abund(x.sample, index.family)
#       results.n[i, ] <- unlist(x.beta)
#       #setTxtProgressBar(pb, i)
#     }
#   }else{
#     stop("function is only implemented for two replicates")
#   }  
#   #close(pb)
#   names(results.n) <- names(x.beta)
#   result <- list(sampled.values = results.n, mean.values = sapply(results.n, 
#                                                                   mean), sd.values = sapply(results.n, sd))
#   return(result)
# }

#plotting the beta-diversity over time
plotPlotTime <- function(indices,measure,#name,
                         x,ylab,xlab,marker) #this function is adjusted to plotting into the Notebook
{
  #pdf(paste0(name,".pdf"),width=17.5/2.54,height=17/2.54,pointsize=7)
  {layout(matrix(1:36,ncol=6,nrow=6,byrow=T))
  par(mar=c(4,4,0.5,0.2),tcl=-0.2,mgp=c(2,0.5,0),lwd=1.5)
  for(i in 1:30) plot(x,indices[(i-1)*5+1:5,measure],type="l",ylim=c(0.1,0.8),col=spcols[i],
                      ylab=ylab,las=1,xlab=xlab,cex.lab=9/7)
  plot(1,type="n",axes=F,ann=F)
  text(1,1,marker,xpd=T,cex=9/7)}
  #dev.off()
}

#dendrograms based on temporal profiles
dendplots <- function(indices,measure,#name,
                      otherdist=env[,c("x","y")],xy=scheme,dcols=spcols) #this function is adjusted to plotting into the Notebook
{
  ppdend <- as.dist(1-cor(unlist(sapply(1:30,function(x) indices[(x-1)*5+1:5,measure]))))
  pdend <- hclust(ppdend,method="ward.D2")
  dend <- as.dendrogram(pdend)
  labels_colors(dend) <- dcols[order.dendrogram(dend)]
  #pdf(paste0(name,".pdf"),width=17.5/2.54,height=6/2.54,pointsize=8)
  {layout(matrix(1:2,ncol=2,nrow=1))
  par(mar=c(1,2.5,0.5,0.5))
  plot(dend)
  mstat <- vector()
  msig <- vector()
  for(i in 0:5){
    m <- mantel(ppdend,
                dist(otherdist[(i*60)+2*1:30,]),"spearman")
    mstat <- append(mstat,m$statistic)
    msig <- append(msig,m$signif)
  }
  par(mar=c(0.5,0.5,0.5,0.5))
  plot(xy$sites,col=dcols,pch=16,cex=5,axes=F,ann=F)
  box()
  ordicluster(xy,pdend,prune=6,col=cutree(pdend,7))
  if(any(msig<0.05)) text(1.5,1.5,round(mean(mstat),digits = 2),col="red") else text(1.5,1.5,round(mean(mstat),digits=2))}
  #dev.off()
  return(list("signif"=msig,"statistic"=mstat))
}

#Mantel-tests to compare dendrograms of temporal profiles to environmental variables (not used in the notebook)
# dendmantel <- function(indices,measure,name,lname="",rname="",otherMeas=env$SoilMois,dcols=spcols,pdf=T)
# {
#   ppdend <- as.dist(1-cor(unlist(sapply(1:30,function(x) indices[(x-1)*5+1:5,measure]))))
#   pdend <- hclust(ppdend,method="ward.D2")
#   dend <- as.dendrogram(pdend)
#   labels_colors(dend) <- dcols[order.dendrogram(dend)]
#   ppsdend <- aggregate(otherMeas,list(rep(1:180,each=2)),mean)[,2]
#   psdend <- as.dist(1-cor(t(unlist(sapply(1:6,function(x) ppsdend[(x-1)*30+1:30])))))
#   sdend <- hclust(psdend,method="ward.D2")
#   dend2 <- as.dendrogram(sdend)
#   labels_colors(dend2) <- dcols[order.dendrogram(dend2)]
#   if(pdf) pdf(paste0(name,".pdf"),width=17.5/2.54,height=6/2.54,pointsize=7)
#   #layout(matrix(1:2,ncol=2,nrow=1))
#   #par(mar=c(1,2.5,0.5,0.5))
#   #plot(dend)
#   mstat <- vector()
#   msig <- vector()
#   for(i in 0:5){
#     m <- mantel(ppdend,
#                 psdend,"spearman")
#     mstat <- append(mstat,m$statistic)
#     msig <- append(msig,m$signif)
#   }
#   par(mar=c(0.2,0.2,0.2,0.2))
#   
#   #plot(xy$sites,col=dcols,pch=16,cex=5,axes=F,ann=F)
#   #box()
#   #ordicluster(xy,pdend,prune=6,col=cutree(pdend,7))
#   if(any(msig<0.05)){
#     tanglegram( dend,  dend2 ,sort=T , sub=round(mean(mstat),digits = 2),main_left=lname,main_right = rname,lwd = 1.5)
#   }else{
#     tanglegram( dend,  dend2 ,sort=T,main_left=lname,main_right = rname ,lwd = 1.5)
#   }
#   if(pdf) dev.off()
#   return(list("signif"=msig,"statistic"=mstat))
# }

#more comparisons of distance matrices (not used in the notebook)
# diff.sample.kezia <- function (x, func = "euc", sites = length(x), tps = 6, repl = 2, samples = 64) 
# {
#   if(samples > repl^tps) {
#     #print(paste("Number of requested samples greater than maximum number of permutations. Max. permutations",
#     #            repl^tps, "used."))
#     samples <- repl^tps
#   }
#   #pb <- txtProgressBar(min = 0, max = samples, style = 3)
#   if (sites > length(x)) 
#     stop("More sites requested for sample than are in the dataset")
#   func <- match.arg(func, c("euc", "cv", "mean"))
#   results.n <- vector("numeric", length(samples))
#   if(repl == 2){
#     sampleVec <- sample(1:(repl^tps),samples)
#     for (i in 1:samples) {
#       pos <- as.integer(intToBits(sampleVec[i]-1))[1:tps]
#       sample.pos <- 2*1:tps-pos
#       x.sample <- x[sample.pos]
#       if(func=="euc"){
#         x.value <- mean(dist(x.sample))
#       }else if(func=="cv"){
#         x.value <- abs(sd(x.sample)/mean(x.sample))
#       }else if(func == "mean"){
#         x.value <- mean(x.sample)
#       }
#       results.n[i] <- unlist(x.value)
#       #setTxtProgressBar(pb, i)
#     }
#   }else{
#     stop("function is only implemented for two replicates")
#   }  
#   #close(pb)
#   #names(results.n) <- names(x.value)
#   result <- list(sampled.values = results.n, mean.values = mean(results.n), sd.values = sd(results.n))
#   return(result)
# }

Analysis of beta-diversity over time

Temporal turnover and beta-diversity

To analyse beta-diversity between timepoints, the Soerensen indices between community structures in a subplot at consecutive time points is measured. In addition to the Soerensen index (“SOR”, total dissimilarity), the Simpson dissimilarity (“SIM”, replacement) and the nestedness fraction of the Soerensen index (“SNE” = SOR-SIM) are reported.

library(betapart)
## Warning: package 'betapart' was built under R version 3.4.3
betaPartIndicesTimeDelta <- matrix(0,nrow=length(unique(subplots))*5,ncol=3,
                                   dimnames=list(paste0(rep(subplots,each=5),"_delta_",1:5),c("beta.SIM","beta.SNE","beta.SOR")))
p <- 1
for(i in subplots){
  for(t in 1:5){
    betaPartIndicesTimeDelta[5*(p-1)+t,] <- unlist(beta.sample.kezia(otuB[tabPlots$subplot==i&tabPlots$tp %in% c(t,t+1),],tps=2)$mean.values)
  }
  p <- p + 1
}

This data was used as basis for Figure 3 and Figure S9. Figure 3 is regenerated here:

library(gplots)
{barplot2(t(as.matrix(aggregate(betaPartIndicesTimeDelta[,3],list(rep(1:5,30)),mean)[,-1])),col="grey80",names.arg = c("Apr-May","May-Jun","Jun-Aug","Aug-Oct","Oct-Nov"),las=1,xlab="beta-diversity",ci.l=t(as.matrix(aggregate(betaPartIndicesTimeDelta[,3],list(rep(1:5,30)),mean)[,-1]))-t(as.matrix(aggregate(betaPartIndicesTimeDelta[,2],list(rep(1:5,30)),sd)[,-1])),ci.u=t(as.matrix(aggregate(betaPartIndicesTimeDelta[,3],list(rep(1:5,30)),mean)[,-1]))+t(as.matrix(aggregate(betaPartIndicesTimeDelta[,2],list(rep(1:5,30)),sd)[,-1])),plot.ci=T,horiz=T)
barplot2(t(as.matrix(aggregate(betaPartIndicesTimeDelta[,1],list(rep(1:5,30)),mean)[,-1])),col=c("grey35"),ci.l=t(as.matrix(aggregate(betaPartIndicesTimeDelta[,1],list(rep(1:5,30)),mean)[,-1]))-t(as.matrix(aggregate(betaPartIndicesTimeDelta[,1],list(rep(1:5,30)),sd)[,-1])),ci.u=t(as.matrix(aggregate(betaPartIndicesTimeDelta[,1],list(rep(1:5,30)),mean)[,-1]))+t(as.matrix(aggregate(betaPartIndicesTimeDelta[,1],list(rep(1:5,30)),sd)[,-1])),plot.ci=T,add=T,ci.color="white",horiz=T,ann=F,axes=F)}

Genus-specific temporal turnover and beta-diversity

The same analysis was performed with subsets of the OTUs that belonged to the most important genera, namely Glomus and Claroidoglomus. In the publication, this data is displayed in Figure S8, regenerated here for Glomus:

betaPartIndicesTimeDelta_glomus <- matrix(0,nrow=length(unique(subplots))*5,ncol=3,
                                   dimnames=list(paste0(rep(subplots,each=5),"_delta_",1:5),c("beta.SIM","beta.SNE","beta.SOR")))
p <- 1
for(i in subplots){
  for(t in 1:5){
    betaPartIndicesTimeDelta_glomus[5*(p-1)+t,] <- unlist(beta.sample.kezia(glomusB[tabPlots$subplot==i&tabPlots$tp %in% c(t,t+1),],tps=2)$mean.values)
  }
  p <- p + 1
}


{barplot2(t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,3],list(rep(1:5,30)),mean)[,-1])),col="lightblue",names.arg = c("Apr-May","May-Jun","Jun-Aug","Aug-Oct","Oct-Nov"),las=1,xlab="beta-diversity",ci.l=t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,3],list(rep(1:5,30)),mean)[,-1]))-t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,2],list(rep(1:5,30)),sd)[,-1])),ci.u=t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,3],list(rep(1:5,30)),mean)[,-1]))+t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,2],list(rep(1:5,30)),sd)[,-1])),plot.ci=T,horiz=T)
barplot2(t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,1],list(rep(1:5,30)),mean)[,-1])),col=c("darkblue"),ci.l=t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,1],list(rep(1:5,30)),mean)[,-1]))-t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,1],list(rep(1:5,30)),sd)[,-1])),ci.u=t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,1],list(rep(1:5,30)),mean)[,-1]))+t(as.matrix(aggregate(betaPartIndicesTimeDelta_glomus[,1],list(rep(1:5,30)),sd)[,-1])),plot.ci=T,add=T,ci.color="white",horiz=T,ann=F,axes=F)}

And for Claroidoglomus:

betaPartIndicesTimeDelta_claro <- matrix(0,nrow=length(unique(subplots))*5,ncol=3,
                                   dimnames=list(paste0(rep(subplots,each=5),"_delta_",1:5),c("beta.SIM","beta.SNE","beta.SOR")))
p <- 1
for(i in subplots){
  for(t in 1:5){
    betaPartIndicesTimeDelta_claro[5*(p-1)+t,] <- unlist(beta.sample.kezia(claroB[tabPlots$subplot==i&tabPlots$tp %in% c(t,t+1),],tps=2)$mean.values)
  }
  p <- p + 1
}


{barplot2(t(as.matrix(aggregate(betaPartIndicesTimeDelta_claro[,3],list(rep(1:5,30)),mean)[,-1])),col="darkolivegreen2",names.arg = c("Apr-May","May-Jun","Jun-Aug","Aug-Oct","Oct-Nov"),las=1,xlab="beta-diversity",ci.l=t(as.matrix(aggregate(betaPartIndicesTimeDelta_claro[,3],list(rep(1:5,30)),mean)[,-1]))-t(as.matrix(aggregate(betaPartIndicesTimeDelta_claro[,2],list(rep(1:5,30)),sd)[,-1])),ci.u=t(as.matrix(aggregate(betaPartIndicesTimeDelta_claro[,3],list(rep(1:5,30)),mean)[,-1]))+t(as.matrix(aggregate(betaPartIndicesTimeDelta_claro[,2],list(rep(1:5,30)),sd)[,-1])),plot.ci=T,horiz=T)
barplot2(t(as.matrix(aggregate(betaPartIndic