#Import libraries and set directory library(vegan) library(MASS) dir <- "M:/OTC" setwd(dir) #import raw OTU table otc_r<-read.table("otu_raw.txt", header=T) head(otc_r) str(otc_r) otc_r<-as.data.frame(t(otc_r)) #import stand table indicating location, treatment, and combination stand<-read.table("stand.txt", header=T) stand$location<-as.factor(stand$location) stand$treatment<-as.factor(stand$treatment) stand$loc_treat<-as.factor(stand$loc_treat) str(stand) #rarefaction curve rare<-min(rowSums(otc_r)) rarecurve(otc_r, step = 20, sample = rare, cex = 0.7, xlab="Reads", ylab="OTUs") #species accumulation curve specacu<-specaccum(otc_r) plot(specacu, ylab="OTUs", xlab="Samples") ##FIGURE 1a and 1b##rarefaction curve and species accumulation curve## #with sample names on rarecurve windows(width=7.00787402, height=3.38582677) par(mfrow=c(1,2),mai=c(0.6,0.75,0.1,0.15), mar=c(3,3,1,1),mgp=c(2,1,0)) plot(specacu, ylab="OTUs", xlab="Samples") legend("topleft", "a", bty="n") rarecurve(otc_r, step = 20, cex = 0.7, xlab="Reads", ylab="OTUs", label=TRUE) legend("topleft", "b", bty="n") #without sample names on rarecurve windows(width=7.00787402, height=3.38582677) par(mfrow=c(1,2),mai=c(0.6,0.75,0.1,0.15), mar=c(3,3,1,1),mgp=c(2,1,0)) plot(specacu, ylab="OTUs", xlab="Samples") legend("topleft", "a", bty="n") rarecurve(otc_r, step = 20, cex = 0.7, xlab="Reads", ylab="OTUs", label=FALSE) legend("topleft", "b", bty="n") ##FIGURE 2a and 2b## GNMDS and CCA of unrarefied pres/abs## #presabs import matrix otc<-read.table("pres_abs.txt", header=T) head(otc) str(otc) otc<-as.data.frame(t(otc)) #run gnmds dist_5<-vegdist(otc,method="bray") mds_5=NULL for(i in 1:100) {mds_5[[i]]<-isoMDS(dist_5,initMDS(dist_5, k=2), k=2, maxit=1000,tol=1e-7)} mds.stress_5<-unlist(lapply(mds_5,function(v){v[[2]]})) ordered_5<-order(mds.stress_5) mds.best_5<-postMDS(mds_5[[ordered_5[1]]],dist_5) mds.secbest_5<-postMDS(mds_5[[ordered_5[2]]],dist_5) procrustes(mds.best_5, mds.secbest_5, permutations=999) protest(mds.best_5, mds.secbest_5, permutations=999) gnmds1_5<-mds.best_5$points[,1] gnmds2_5<-mds.best_5$points[,2] #import stand stable stand<-read.table("stand.txt", header=T) stand$location<-as.factor(stand$location) stand$treatment<-as.factor(stand$treatment) stand$loc_treat<-as.factor(stand$loc_treat) str(stand) #import environmental variables table env<-read.table("env.txt", header=T) str(env) env$location<-as.factor(env$location) env$treatment<-as.factor(env$treatment) str(env) envfit_5<-envfit(gnmds1_5~env$location+env$treatment+env$richness+env$root_weight, permutations=999) envfit_5 #CCA presence/absence cca<-cca(otc~env$treatment+Condition(env$location)) cca testcca<-permutest(cca,permutations=999) testcca testcca_a<-anova.cca(cca,permutations=999) testcca_a #variation pertitioning varp <- varpart (otc, ~ treatment, ~ location, ~ root_weight, ~ richness, data = env) varp plot (varp, digits = 2) #Fig 2a and 2b plot windows(width=6.85039, height=3.3) par(mfrow=c(1,2),mai=c(0.6,0.75,0.1,0.15), mar=c(2.8,2.8,1,1),mgp=c(2,1,0),oma = c(1, 1, 1, 4)) #Fig 2a GNMDS unrarefied pres/abs plot(gnmds1_5, gnmds2_5, xlab="GNMDS axis 1", ylab="GNMDS axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(gnmds1_5[stand$loc_treat==1],gnmds2_5[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(gnmds1_5[stand$loc_treat==2],gnmds2_5[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(gnmds1_5[stand$loc_treat==3],gnmds2_5[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(gnmds1_5[stand$loc_treat==4],gnmds2_5[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(gnmds1_5[stand$loc_treat==5],gnmds2_5[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(gnmds1_5[stand$loc_treat==6],gnmds2_5[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("topleft", "a", bty="n") pchvec <- c(16,16, 16,17, 17,17) colvec <- c("#1b9e77","#d95f02","#7570b3","#1b9e77","#d95f02","#7570b3") #Fig 2b CCA unrarefied pres/abs plot(cca, xlab="CCA axis 1", ylab="CCA axis 2",type="n", cex.axis=0.8, cex.lab=0.9, choices=c(1,2),display="sites") with(stand, points(cca, display = "sites", pch=pchvec[loc_treat], col=colvec[loc_treat], cex=2)) legend("topleft","b", bty="n") par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE) #legend plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n") legend("right", c("OTC 1","Control 1","OTC 2","Control 2","OTC 3","Control 3"), xpd = TRUE, horiz = FALSE, inset = c(0, 0), bty = "n", pch=c(16,17,16,17,16,17), col=c("#1b9e77","#1b9e77","#d95f02","#d95f02","#7570b3","#7570b3"), cex = 0.7) ##FIGURE 3a and 3b##order bar plot and genera within Agaricales barplot## #Order table bar plot order<-read.table("order5.txt", header=T) head(order) str(order) sumsorder<-colSums(order) T1<-order[,1]/sumsorder[1] T2<-order[,2]/sumsorder[2] T3<-order[,3]/sumsorder[3] C1<-order[,4]/sumsorder[4] C2<-order[,5]/sumsorder[5] T4<-order[,6]/sumsorder[6] T5<-order[,7]/sumsorder[7] T6<-order[,8]/sumsorder[8] C3<-order[,9]/sumsorder[9] C4<-order[,10]/sumsorder[10] T7<-order[,11]/sumsorder[11] T8<-order[,12]/sumsorder[12] T9<-order[,13]/sumsorder[13] C5<-order[,14]/sumsorder[14] C6<-order[,15]/sumsorder[15] orderper<-cbind(T1,T2,T3,C1,C2,T4,T5,T6,C3,C4,T7,T8,T9,C5,C6) rownames(orderper)<-rownames(order) colnames(orderper)<-colnames(order) orderper<-as.matrix(orderper) orderper barplot(orderper) windows(height = 3.38582,width = 6.85039) par(mgp=c(2,1,0), mar=c(3.5,3.5,2,7),fg="black") barplot(orderper,col=c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a","#ffff99","#b15928","#13496d","darksalmon"),cex.axis=0.8, cex.names=0.8) par(fg="white") legend("topright",inset = c(-0.28,0), legend=rownames(order),xpd=TRUE, pch=0,bty="n",cex=0.8, col=1) par(fg="black") legend("topright",inset = c(-0.28,0), legend=rownames(order),xpd=TRUE, pch=15,bty="n",cex=0.8, col=c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a","#ffff99","#b15928","#13496d","darksalmon")) #agaricales table bar plot agaric<-read.table("agaric.txt", header=T) head(agaric) str(agaric) sumsagaric<-colSums(agaric) T1<-agaric[,1]/sumsagaric[1] T2<-agaric[,2]/sumsagaric[2] T3<-agaric[,3]/sumsagaric[3] C1<-agaric[,4]/sumsagaric[4] C2<-agaric[,5]/sumsagaric[5] T4<-agaric[,6]/sumsagaric[6] T5<-agaric[,7]/sumsagaric[7] T6<-agaric[,8]/sumsagaric[8] C3<-agaric[,9]/sumsagaric[9] C4<-agaric[,10]/sumsagaric[10] T7<-agaric[,11]/sumsagaric[11] T8<-agaric[,12]/sumsagaric[12] T9<-agaric[,13]/sumsagaric[13] C5<-agaric[,14]/sumsagaric[14] C6<-agaric[,15]/sumsagaric[15] agaricper<-cbind(T1,T2,T3,C1,C2,T4,T5,T6,C3,C4,T7,T8,T9,C5,C6) rownames(agaricper)<-rownames(agaric) colnames(agaricper)<-colnames(agaric) agaricper<-as.matrix(agaricper) agaricper barplot(agaricper) windows(height = 3.38582,width = 6.85039) par(mgp=c(2,1,0), mar=c(3.5,3.5,2,7),fg="black") barplot(agaricper,col=c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6"),cex.axis=0.8, cex.names=0.8) par(fg="white") legend("topright",inset = c(-0.28,0), legend=rownames(agaric),xpd=TRUE, pch=0,bty="n",cex=0.8, col=1) par(fg="black") legend("topright",inset = c(-0.28,0), legend=rownames(agaric),xpd=TRUE, pch=15,bty="n",cex=0.8, col=c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6")) #Fig. 3a and b ##windows(width=6.85039, height=3.3) ##par(mfrow=c(1,2),mai=c(0.6,0.75,0.1,0.15), mar=c(2.8,2.8,1,1),mgp=c(2,1,0),oma = c(1, 1, 1, 4)) windows(height = 6.77164,width = 6.85039) par(mfrow=c(2,1),mgp=c(2,1,0), mar=c(3.5,3.5,2,7),fg="black") barplot(orderper,col=c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a","#ffff99","#b15928","#13496d","darksalmon"),cex.axis=0.8, cex.names=0.8) par(fg="white") legend("topright",inset = c(-0.28,0), legend=rownames(order),xpd=TRUE, pch=0,bty="n",cex=0.8, col=1) par(fg="black") legend("topright",inset = c(-0.28,0), legend=rownames(order),xpd=TRUE, pch=15,bty="n",cex=0.8, col=c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6","#6a3d9a","#ffff99","#b15928","#13496d","darksalmon")) barplot(agaricper,col=c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6"),cex.axis=0.8, cex.names=0.8) par(fg="white") legend("topright",inset = c(-0.22,0), legend=rownames(agaric),xpd=TRUE, pch=0,bty="n",cex=0.8, col=1) par(fg="black") legend("topright",inset = c(-0.22,0), legend=rownames(agaric),xpd=TRUE, pch=15,bty="n",cex=0.8, col=c("#a6cee3","#1f78b4","#b2df8a","#33a02c","#fb9a99","#e31a1c","#fdbf6f","#ff7f00","#cab2d6")) par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE) plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n") legend("topleft","a", bty="n") par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE) plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n") legend("left","b", bty="n") ####SUPPLEMENT otc<-read.table("pres_abs.txt", header=T) head(otc) str(otc) otc<-as.data.frame(t(otc)) dist_5<-vegdist(otc,method="bray") mds_5=NULL for(i in 1:100) {mds_5[[i]]<-isoMDS(dist_5,initMDS(dist_5, k=2), k=2, maxit=1000,tol=1e-7)} mds.stress_5<-unlist(lapply(mds_5,function(v){v[[2]]})) ordered_5<-order(mds.stress_5) mds.best_5<-postMDS(mds_5[[ordered_5[1]]],dist_5) mds.secbest_5<-postMDS(mds_5[[ordered_5[2]]],dist_5) procrustes(mds.best_5, mds.secbest_5, permutations=999) protest(mds.best_5, mds.secbest_5, permutations=999) gnmds1_5<-mds.best_5$points[,1] gnmds2_5<-mds.best_5$points[,2] dca_5<- decorana(otc) dca1_5<-scores(dca_5, display="sites", origin=FALSE)[,1] dca2_5<-scores(dca_5, display="sites", origin=FALSE)[,2] dca3_5<-scores(dca_5, display="sites", origin=FALSE)[,3] cor.test(gnmds1_5,dca2_5,method="kendall") cor.test(gnmds1_5,dca1_5,method="kendall") cor.test(gnmds2_5,dca2_5,method="kendall") stand<-read.table("stand.txt", header=T) stand$location<-as.factor(stand$location) stand$treatment<-as.factor(stand$treatment) stand$loc_treat<-as.factor(stand$loc_treat) str(stand) ########DCA and GNMDS plots Pres/Abs plot(dca1_5, dca2_5, xlab="DCA axis 1", ylab="DCA axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(dca1_5[stand$loc_treat==1],dca2_5[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(dca1_5[stand$loc_treat==2],dca2_5[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(dca1_5[stand$loc_treat==3],dca2_5[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(dca1_5[stand$loc_treat==4],dca2_5[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(dca1_5[stand$loc_treat==5],dca2_5[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(dca1_5[stand$loc_treat==6],dca2_5[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("topleft",c("OTC 1","Control OTC 1","OTC 2","Control OTC 2","OTC 3","Control OTC 3"), col=c("#1b9e77","#1b9e77","#d95f02","#d95f02","#7570b3","#7570b3"), pch=c(16,17,16,17,16,17), cex=1.5) plot(gnmds1_5, gnmds2_5, xlab="GNMDS axis 1", ylab="GNMDS axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(gnmds1_5[stand$loc_treat==1],gnmds2_5[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(gnmds1_5[stand$loc_treat==2],gnmds2_5[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(gnmds1_5[stand$loc_treat==3],gnmds2_5[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(gnmds1_5[stand$loc_treat==4],gnmds2_5[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(gnmds1_5[stand$loc_treat==5],gnmds2_5[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(gnmds1_5[stand$loc_treat==6],gnmds2_5[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("topright",c("OTC 1","Control OTC 1","OTC 2","Control OTC 2","OTC 3","Control OTC 3"), col=c("#1b9e77","#1b9e77","#d95f02","#d95f02","#7570b3","#7570b3"), pch=c(16,17,16,17,16,17), cex=1.0) #### #Raw Reads otc_r<-read.table("otu_raw.txt", header=T) head(otc_r) str(otc_r) otc_r<-as.data.frame(t(otc_r)) dist_r<-vegdist(otc_r,method="bray") mds_r=NULL for(i in 1:100) {mds_r[[i]]<-isoMDS(dist_r,initMDS(dist_r, k=2), k=2, maxit=1000,tol=1e-7)} mds.stress_r<-unlist(lapply(mds_r,function(v){v[[2]]})) ordered_r<-order(mds.stress_r) mds.best_r<-postMDS(mds_r[[ordered_r[1]]],dist_r) mds.secbest_r<-postMDS(mds_r[[ordered_r[2]]],dist_r) procrustes(mds.best_r, mds.secbest_r, permutations=999) protest(mds.best_r, mds.secbest_r, permutations=999) gnmds1_r<-mds.best_r$points[,1] gnmds2_r<-mds.best_r$points[,2] dca_r<- decorana(otc_r) dca1_r<-scores(dca_r, display="sites", origin=FALSE)[,1] dca2_r<-scores(dca_r, display="sites", origin=FALSE)[,2] dca3_r<-scores(dca_r, display="sites", origin=FALSE)[,3] cor.test(gnmds1_r,dca2_r,method="kendall") cor.test(gnmds1_r,dca1_r,method="kendall") cor.test(gnmds2_r,dca2_r,method="kendall") stand<-read.table("stand.txt", header=T) stand$location<-as.factor(stand$location) stand$treatment<-as.factor(stand$treatment) stand$loc_treat<-as.factor(stand$loc_treat) str(stand) ######Raw Reads Plots plot(gnmds1_r, gnmds2_r, xlab="GNMDS axis 1", ylab="GNMDS axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(gnmds1_r[stand$loc_treat==1],gnmds2_r[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(gnmds1_r[stand$loc_treat==2],gnmds2_r[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(gnmds1_r[stand$loc_treat==3],gnmds2_r[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(gnmds1_r[stand$loc_treat==4],gnmds2_r[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(gnmds1_r[stand$loc_treat==5],gnmds2_r[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(gnmds1_r[stand$loc_treat==6],gnmds2_r[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("topright",c("OTC 1","Control OTC 1","OTC 2","Control OTC 2","OTC 3","Control OTC 3"), col=c("#1b9e77","#1b9e77","#d95f02","#d95f02","#7570b3","#7570b3"), pch=c(16,17,16,17,16,17), cex=1.0) plot(dca1_r, dca2_r, xlab="DCA axis 1", ylab="DCA axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(dca1_r[stand$loc_treat==1],dca2_r[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(dca1_r[stand$loc_treat==2],dca2_r[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(dca1_r[stand$loc_treat==3],dca2_r[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(dca1_r[stand$loc_treat==4],dca2_r[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(dca1_r[stand$loc_treat==5],dca2_r[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(dca1_r[stand$loc_treat==6],dca2_r[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("topleft",c("OTC 1","Control OTC 1","OTC 2","Control OTC 2","OTC 3","Control OTC 3"), col=c("#1b9e77","#1b9e77","#d95f02","#d95f02","#7570b3","#7570b3"), pch=c(16,17,16,17,16,17), cex=1.5, bty="n") ## Rarefied to 1059 otc_rare<-read.table("otu_rare.txt", header=T) head(otc_rare) str(otc_rare) otc_rare<-as.data.frame(t(otc_rare)) dist_sm<-vegdist(otc_rare,method="bray") mds_sm=NULL for(i in 1:100) {mds_sm[[i]]<-isoMDS(dist_sm,initMDS(dist_sm, k=2), k=2, maxit=1000,tol=1e-7)} mds.stress_sm<-unlist(lapply(mds_sm,function(v){v[[2]]})) ordered_sm<-order(mds.stress_sm) mds.best_sm<-postMDS(mds_sm[[ordered_sm[1]]],dist_sm) mds.secbest_sm<-postMDS(mds_sm[[ordered_sm[2]]],dist_sm) procrustes(mds.best_sm, mds.secbest_sm, permutations=999) protest(mds.best_sm, mds.secbest_sm, permutations=999) gnmds1_sm<-mds.best_sm$points[,1] gnmds2_sm<-mds.best_sm$points[,2] dca_sm<- decorana(otc_rare) dca1_sm<-scores(dca_sm, display="sites", origin=FALSE)[,1] dca2_sm<-scores(dca_sm, display="sites", origin=FALSE)[,2] dca3_sm<-scores(dca_sm, display="sites", origin=FALSE)[,3] cor.test(gnmds1_sm,dca2_sm,method="kendall") cor.test(gnmds1_sm,dca1_sm,method="kendall") cor.test(gnmds2_sm,dca2_sm,method="kendall") plot(gnmds1_sm, gnmds2_sm, xlab="GNMDS axis 1", ylab="GNMDS axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(gnmds1_sm[stand$loc_treat==1],gnmds2_sm[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(gnmds1_sm[stand$loc_treat==2],gnmds2_sm[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(gnmds1_sm[stand$loc_treat==3],gnmds2_sm[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(gnmds1_sm[stand$loc_treat==4],gnmds2_sm[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(gnmds1_sm[stand$loc_treat==5],gnmds2_sm[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(gnmds1_sm[stand$loc_treat==6],gnmds2_sm[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("top",c("OTC 1","Control OTC 1","OTC 2","Control OTC 2","OTC 3","Control OTC 3"), col=c("#1b9e77","#1b9e77","#d95f02","#d95f02","#7570b3","#7570b3"), pch=c(16,17,16,17,16,17), cex=1.0, bty="n") plot(dca1_sm, dca2_sm, xlab="DCA axis 1", ylab="DCA axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(dca1_sm[stand$loc_treat==1],dca2_sm[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(dca1_sm[stand$loc_treat==2],dca2_sm[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(dca1_sm[stand$loc_treat==3],dca2_sm[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(dca1_sm[stand$loc_treat==4],dca2_sm[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(dca1_sm[stand$loc_treat==5],dca2_sm[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(dca1_sm[stand$loc_treat==6],dca2_sm[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("topleft",c("OTC 1","Control OTC 1","OTC 2","Control OTC 2","OTC 3","Control OTC 3"), col=c("#1b9e77","#1b9e77","#d95f02","#d95f02","#7570b3","#7570b3"), pch=c(16,17,16,17,16,17), cex=1.0) ####plots together windows(height = 20.31492,width = 13.70079) par(mfrow=c(3,2), mgp=c(2,1,0),mar=c(5,5,3,3)) #raw reads plot(gnmds1_r, gnmds2_r, xlab="GNMDS axis 1", ylab="GNMDS axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(gnmds1_r[stand$loc_treat==1],gnmds2_r[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(gnmds1_r[stand$loc_treat==2],gnmds2_r[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(gnmds1_r[stand$loc_treat==3],gnmds2_r[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(gnmds1_r[stand$loc_treat==4],gnmds2_r[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(gnmds1_r[stand$loc_treat==5],gnmds2_r[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(gnmds1_r[stand$loc_treat==6],gnmds2_r[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("bottomleft","a", bty="n", cex=1.5) plot(dca1_r, dca2_r, xlab="DCA axis 1", ylab="DCA axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(dca1_r[stand$loc_treat==1],dca2_r[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(dca1_r[stand$loc_treat==2],dca2_r[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(dca1_r[stand$loc_treat==3],dca2_r[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(dca1_r[stand$loc_treat==4],dca2_r[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(dca1_r[stand$loc_treat==5],dca2_r[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(dca1_r[stand$loc_treat==6],dca2_r[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("bottomleft","b", bty="n", cex=1.5) #presabs plot(gnmds1_5, gnmds2_5, xlab="GNMDS axis 1", ylab="GNMDS axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(gnmds1_5[stand$loc_treat==1],gnmds2_5[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(gnmds1_5[stand$loc_treat==2],gnmds2_5[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(gnmds1_5[stand$loc_treat==3],gnmds2_5[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(gnmds1_5[stand$loc_treat==4],gnmds2_5[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(gnmds1_5[stand$loc_treat==5],gnmds2_5[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(gnmds1_5[stand$loc_treat==6],gnmds2_5[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("bottomleft","c", bty="n", cex=1.5) plot(dca1_5, dca2_5, xlab="DCA axis 1", ylab="DCA axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(dca1_5[stand$loc_treat==1],dca2_5[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(dca1_5[stand$loc_treat==2],dca2_5[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(dca1_5[stand$loc_treat==3],dca2_5[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(dca1_5[stand$loc_treat==4],dca2_5[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(dca1_5[stand$loc_treat==5],dca2_5[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(dca1_5[stand$loc_treat==6],dca2_5[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("topleft",c("OTC 1","Control OTC 1","OTC 2","Control OTC 2","OTC 3","Control OTC 3"), col=c("#1b9e77","#1b9e77","#d95f02","#d95f02","#7570b3","#7570b3"), pch=c(16,17,16,17,16,17), cex=1.2, bty="n") legend("bottomleft","d", bty="n", cex=1.5) #rarefied plot(gnmds1_sm, gnmds2_sm, xlab="GNMDS axis 1", ylab="GNMDS axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(gnmds1_sm[stand$loc_treat==1],gnmds2_sm[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(gnmds1_sm[stand$loc_treat==2],gnmds2_sm[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(gnmds1_sm[stand$loc_treat==3],gnmds2_sm[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(gnmds1_sm[stand$loc_treat==4],gnmds2_sm[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(gnmds1_sm[stand$loc_treat==5],gnmds2_sm[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(gnmds1_sm[stand$loc_treat==6],gnmds2_sm[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("bottomleft","e", bty="n", cex=1.5) plot(dca1_sm, dca2_sm, xlab="DCA axis 1", ylab="DCA axis 2",type="n", cex.axis=0.8, cex.lab=0.9) points(dca1_sm[stand$loc_treat==1],dca2_sm[stand$loc_treat==1], pch=16, col="#1b9e77", cex=2) points(dca1_sm[stand$loc_treat==2],dca2_sm[stand$loc_treat==2], pch=16, col="#d95f02", cex=2) points(dca1_sm[stand$loc_treat==3],dca2_sm[stand$loc_treat==3], pch=16, col="#7570b3", cex=2) points(dca1_sm[stand$loc_treat==4],dca2_sm[stand$loc_treat==4], pch=17, col="#1b9e77", cex=2) points(dca1_sm[stand$loc_treat==5],dca2_sm[stand$loc_treat==5], pch=17, col="#d95f02", cex=2) points(dca1_sm[stand$loc_treat==6],dca2_sm[stand$loc_treat==6], pch=17, col="#7570b3", cex=2) legend("bottomleft","f", bty="n", cex=1.5)