library(AER) library(multcompView) library(multcomp) library(ggpubr) library(Rmisc) library(ggplot2) library(do) library(grid) # ANOVA in Arable land (A) ------------------------------------------------------------------- library(car) recoll<- read.csv("Coll_recov_all180_nz.csv", header = TRUE) recoll$Treatment = as.factor(recoll$Treatment) recoll$Time= as.factor(recoll$Time) recoll$Time <- ordered(recoll$Time,levels = c("two weeks", "seven weeks", "twelve weeks")) recoll$Total = rowSums(recoll[,6:19]) recoll$lgtotal <- log10(recoll$Total+1) #subset Arable land (A) HaT1<-subset(recoll, Time=="two weeks" & Destination=="A") HaT2<-subset(recoll, Time=="seven weeks" & Destination=="A") HaT3<-subset(recoll, Time=="twelve weeks" & Destination=="A") #shapiro.test of lgtogal, when p>0.05, sample is not normally distributed shapiro.test(HaT1[1:5,21]) shapiro.test(HaT1[6:10,21]) shapiro.test(HaT1[11:15,21]) shapiro.test(HaT1[16:20,21]) shapiro.test(HaT2[1:5,21]) shapiro.test(HaT2[6:10,21]) shapiro.test(HaT2[11:15,21]) shapiro.test(HaT2[16:20,21]) shapiro.test(HaT3[1:5,21]) shapiro.test(HaT3[6:10,21]) shapiro.test(HaT3[11:15,21]) shapiro.test(HaT3[16:20,21]) #leveneTest leveneTest(lgtotal ~ Treatment, data=HaT1) leveneTest(lgtotal ~ Treatment, data=HaT2) leveneTest(lgtotal ~ Treatment, data=HaT3) # Anova, lgtotal ~ Treatment, at three time points in "A" aov.ha1 <- aov(lgtotal ~ Treatment, data = HaT1) summary(aov.ha1) TukeyHSD(aov.ha1) aov.ha2 <- aov(lgtotal ~ Treatment, data = HaT2) summary(aov.ha2) TukeyHSD(aov.ha2) aov.ha3 <- aov(lgtotal ~ Treatment, data = HaT3) summary(aov.ha3) TukeyHSD(aov.ha3) # Anova, species number ~ Treatment, at three time points in "A" aov.s.ha1 <- aov(Spnumber ~ Treatment, data = HaT1) summary(aov.s.ha1) TukeyHSD(aov.s.ha1) aov.s.ha2 <- aov(Spnumber ~ Treatment, data = HaT2) summary(aov.s.ha2) TukeyHSD(aov.s.ha2) aov.s.ha3 <- aov(Spnumber ~ Treatment, data = HaT3) summary(aov.s.ha3) TukeyHSD(aov.s.ha3) # ANOVA in Grassland (G) -------------------------------------------------------------- #subset habitats G HgT1<-subset(recoll, Time=="two weeks" & Destination=="G") HgT2<-subset(recoll, Time=="seven weeks" & Destination=="G") HgT3<-subset(recoll, Time=="twelve weeks" & Destination=="G") #shapiro.test of lgtogal, when p>0.05, sample is not normally distributed shapiro.test(HgT1[1:5,21]) shapiro.test(HgT1[6:10,21]) shapiro.test(HgT1[11:15,21]) shapiro.test(HgT1[16:20,21]) shapiro.test(HgT2[1:5,21]) shapiro.test(HgT2[6:10,21]) shapiro.test(HgT2[11:15,21]) shapiro.test(HgT2[16:20,21]) shapiro.test(HgT3[1:5,21]) shapiro.test(HgT3[6:10,21]) shapiro.test(HgT3[11:15,21]) shapiro.test(HgT3[16:20,21]) #leveneTest leveneTest(lgtotal ~ Treatment, data=HgT1) leveneTest(lgtotal ~ Treatment, data=HgT2) leveneTest(lgtotal ~ Treatment, data=HgT3) # Anova, lgtotal ~ Treatment, at three time points in "G" aov.hg1 <- aov(lgtotal ~ Treatment, data = HgT1) summary(aov.hg1) TukeyHSD(aov.hg1) aov.hg2 <- aov(lgtotal ~ Treatment, data = HgT2) summary(aov.hg2) TukeyHSD(aov.hg2) aov.hg3 <- aov(lgtotal ~ Treatment, data = HgT3) summary(aov.hg3) TukeyHSD(aov.hg3) # Anova, species number ~ Treatment, at three time points in "G" aov.s.hg1 <- aov(Spnumber ~ Treatment, data = HgT1) summary(aov.s.hg1) TukeyHSD(aov.s.hg1) aov.s.hg2 <- aov(Spnumber ~ Treatment, data = HgT2) summary(aov.s.hg2) TukeyHSD(aov.s.hg2) aov.s.hg3 <- aov(Spnumber ~ Treatment, data = HgT3) summary(aov.s.hg3) TukeyHSD(aov.s.hg3) # ANOVA in BSAL -------------------------------------------------------------- #subset habitats B HbT1<-subset(recoll, Time=="two weeks" & Destination=="B") HbT2<-subset(recoll, Time=="seven weeks" & Destination=="B") HbT3<-subset(recoll, Time=="twelve weeks" & Destination=="B") #shapiro.test of lgtogal, when p>0.05, sample is not normally distributed shapiro.test(HbT1[1:5,21]) shapiro.test(HbT1[6:10,21]) shapiro.test(HbT1[11:15,21]) shapiro.test(HbT1[16:20,21]) shapiro.test(HbT2[1:5,21]) shapiro.test(HbT2[6:10,21]) shapiro.test(HbT2[11:15,21]) shapiro.test(HbT2[16:20,21]) shapiro.test(HbT3[1:5,21]) shapiro.test(HbT3[6:10,21]) shapiro.test(HbT3[11:15,21]) shapiro.test(HbT3[16:20,21]) #leveneTest leveneTest(lgtotal ~ Treatment, data=HbT1) leveneTest(lgtotal ~ Treatment, data=HbT2) leveneTest(lgtotal ~ Treatment, data=HbT3) # Anova, lgtotal ~ Treatment, at three time points in "G" aov.hb1 <- aov(lgtotal ~ Treatment, data = HbT1) summary(aov.hb1) TukeyHSD(aov.hb1) aov.hb2 <- aov(lgtotal ~ Treatment, data = HbT2) summary(aov.hb2) TukeyHSD(aov.hb2) aov.hb3 <- aov(lgtotal ~ Treatment, data = HbT3) summary(aov.hb3) TukeyHSD(aov.hb3) # Anova, species number ~ Treatment, at three time points in "G" aov.s.hb1 <- aov(Spnumber ~ Treatment, data = HbT1) summary(aov.s.hb1) TukeyHSD(aov.s.hb1) aov.s.hb2 <- aov(Spnumber ~ Treatment, data = HbT2) summary(aov.s.hb2) TukeyHSD(aov.s.hb2) aov.s.hb3 <- aov(Spnumber ~ Treatment, data = HbT3) summary(aov.s.hb3) TukeyHSD(aov.s.hb3) # Figures abundance ----------------------------------- # in Arable land recoll<- read.csv("Coll_recov_all180_nz.csv", header = TRUE) recoll$Treatment = as.factor(recoll$Treatment) recoll$Time= as.factor(recoll$Time) recoll$Time <- ordered(recoll$Time,levels = c("two weeks", "seven weeks", "twelve weeks")) recoll$Total = rowSums(recoll[,6:19]) recoll<-tidyr::unite(recoll, "site", Time, Treatment, Destination, sep = "",remove=FALSE) recoll$lgtotal <- log10(recoll$Total+1) Ha<-subset(recoll, Destination=="A") Hg<-subset(recoll, Destination=="G") Hb<-subset(recoll, Destination=="B") summary(Ha) Haabu <- summarySE(Ha, measurevar="lgtotal", groupvars=c("Time","Treatment")) Haabu$Time = gsub("two weeks", "2", Haabu$Time) Haabu$Time = gsub("seven weeks", "7", Haabu$Time) Haabu$Time = gsub("twelve weeks","12", Haabu$Time) Haabu$Time = as.factor(Haabu$Time) Haabu$Time = ordered(Haabu$Time,levels = c("2", "7", "12")) Haabu$Treatment= as.factor(Haabu$Treatment) Haabu$Treatment <- ordered(Haabu$Treatment,levels = c("WAA", "OAA", "OGA", "OBA")) summary(Haabu) levels(Haabu$Time) label_Haabu <- c("", "", "", "", "c", "c", "a", "b", "", "","", "") Haabu_fig <- ggplot(Haabu, aes(x=Time, y=lgtotal,color=Treatment, fill=Treatment)) + geom_bar(position=position_dodge(), stat="identity")+ geom_errorbar(aes(ymin=lgtotal-se, ymax=lgtotal+se), colour="black", width=.1, position=position_dodge(.9)) + geom_text(aes(y = lgtotal + 1.8 * se, label = label_Haabu, group = Treatment), colour="red", position = position_dodge(0.8), size = 3, fontface = "bold") + ggtitle("A") + xlab("Sampling time (weeks)") + ylab(" lg(total abundance per block +1)") + scale_y_continuous(expand = c(0, 0)) + geom_text(aes(y = 3), label = "") + theme_bw() Haabu_fig # in Grassland (G) Hgabu <- summarySE(Hg, measurevar="lgtotal", groupvars=c("Time","Treatment")) Hgabu$Time = gsub("two weeks", "2", Hgabu$Time) Hgabu$Time = gsub("seven weeks", "7", Hgabu$Time) Hgabu$Time = gsub("twelve weeks","12", Hgabu$Time) Hgabu$Time = as.factor(Hgabu$Time) Hgabu$Time = ordered(Hgabu$Time,levels = c("2", "7", "12")) Hgabu$Treatment= as.factor(Hgabu$Treatment) Hgabu$Treatment <- ordered(Hgabu$Treatment,levels = c("WGG", "OAG", "OGG", "OBG")) summary(Hgabu) levels(Hgabu$Time) label_Hgabu <- c("ab","b","a","ab","b","ab","a","ab","","","","") Hgabu_fig <- ggplot(Hgabu, aes(x=Time, y=lgtotal,color=Treatment, fill=Treatment)) + geom_bar(position=position_dodge(), stat="identity")+ geom_errorbar(aes(ymin=lgtotal-se, ymax=lgtotal+se), colour="black", width=.1, position=position_dodge(.9)) + geom_text(aes(y = lgtotal + 1.8 * se, label = label_Hgabu, group = Treatment), colour="red", position = position_dodge(0.8), size = 3, fontface = "bold") + ggtitle("B") + xlab("Sampling time (weeks)") + ylab("lg(total abundance per block +1)") + scale_y_continuous(expand = c(0, 0)) + geom_text(aes(y = 3), label = "") + theme_bw() Hgabu_fig # in BSAL (B) Hbabu <- summarySE(Hb, measurevar="lgtotal", groupvars=c("Time","Treatment")) Hbabu$Time = gsub("two weeks", "2", Hbabu$Time) Hbabu$Time = gsub("seven weeks", "7", Hbabu$Time) Hbabu$Time = gsub("twelve weeks","12", Hbabu$Time) Hbabu$Time = as.factor(Hbabu$Time) Hbabu$Time = ordered(Hbabu$Time,levels = c("2", "7", "12")) Hbabu$Treatment= as.factor(Hbabu$Treatment) Hbabu$Treatment <- ordered(Hbabu$Treatment,levels = c("WBB", "OAB", "OGB", "OBB")) summary(Hbabu) levels(Hbabu$Time) label_Hbabu <- c("ab", "b", "a", "b", "b", "c", "a", "b", "", "","", "") Hbabu_fig <- ggplot(Hbabu, aes(x=Time, y=lgtotal,color=Treatment, fill=Treatment)) + geom_bar(position=position_dodge(), stat="identity")+ geom_errorbar(aes(ymin=lgtotal-se, ymax=lgtotal+se), colour="black", width=.1, position=position_dodge(.9)) + geom_text(aes(y = lgtotal + 1.8 * se, label = label_Hbabu, group = Treatment), colour="red", position = position_dodge(0.8), size = 3, fontface = "bold") + ggtitle("C") + xlab("Sampling time (weeks)") + ylab("lg(total abundance per block +1)") + scale_y_continuous(expand = c(0, 0)) + geom_text(aes(y = 3), label = "") + theme_bw() Hbabu_fig grid.newpage() pushViewport(viewport(layout = grid.layout(1, 3))) vplayout = function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) print(Haabu_fig, vp = vplayout(1, 1)) print(Hgabu_fig, vp = vplayout(1, 2)) print(Hbabu_fig, vp = vplayout(1, 3)) # Figures species --------------------------------------------------------- #in Arable land (A) Hasp <- summarySE(Ha, measurevar="Spnumber", groupvars=c("Time","Treatment")) Hasp$Time = gsub("two weeks", "2", Hasp$Time) Hasp$Time = gsub("seven weeks", "7", Hasp$Time) Hasp$Time = gsub("twelve weeks","12", Hasp$Time) Hasp$Time = as.factor(Hasp$Time) Hasp$Time = ordered(Hasp$Time,levels = c("2", "7", "12")) Hasp$Treatment= as.factor(Hasp$Treatment) Hasp$Treatment <- ordered(Hasp$Treatment,levels = c("WAA", "OAA", "OGA", "OBA")) summary(Hasp) levels(Hasp$Time) label_Hasp <- c("", "", "", "", "", "", "", "", "", "","", "") Hasp_fig <- ggplot(Hasp, aes(x=Time, y=Spnumber,color=Treatment, fill=Treatment)) + geom_bar(position=position_dodge(), stat="identity")+ geom_errorbar(aes(ymin=Spnumber-se, ymax=Spnumber+se), colour="black", width=.1, position=position_dodge(.9)) + geom_text(aes(y = Spnumber + 1.8 * se, label = label_Hasp, group = Treatment), colour="red", position = position_dodge(0.8), size = 3, fontface = "bold") + ggtitle("A") + xlab("Sampling time (weeks)") + ylab(" species number per soil block") + scale_y_continuous(expand = c(0, 0)) + geom_text(aes(y = 12), label = "") + theme_bw() Hasp_fig #in Grassland (G) Hgsp <- summarySE(Hg, measurevar="Spnumber", groupvars=c("Time","Treatment")) Hgsp$Time = gsub("two weeks", "2", Hgsp$Time) Hgsp$Time = gsub("seven weeks", "7", Hgsp$Time) Hgsp$Time = gsub("twelve weeks","12", Hgsp$Time) Hgsp$Time = as.factor(Hgsp$Time) Hgsp$Time = ordered(Hgsp$Time,levels = c("2", "7", "12")) Hgsp$Treatment= as.factor(Hgsp$Treatment) Hgsp$Treatment <- ordered(Hgsp$Treatment,levels = c("WGG", "OAG", "OGG", "OBG")) summary(Hgsp) levels(Hgsp$Time) label_Hgsp <- c("", "", "", "", "", "", "", "", "", "","", "") Hgsp_fig <- ggplot(Hgsp, aes(x=Time, y=Spnumber,color=Treatment, fill=Treatment)) + geom_bar(position=position_dodge(), stat="identity")+ geom_errorbar(aes(ymin=Spnumber-se, ymax=Spnumber+se), colour="black", width=.1, position=position_dodge(.9)) + geom_text(aes(y = Spnumber + 1.8 * se, label = label_Hgsp, group = Treatment), colour="red", position = position_dodge(0.8), size = 3, fontface = "bold") + ggtitle("B") + xlab("Sampling time (weeks)") + ylab(" species number per soil block") + scale_y_continuous(expand = c(0, 0)) + geom_text(aes(y = 12), label = "") + theme_bw() Hgsp_fig #in BSAL (B) Hbsp <- summarySE(Hb, measurevar="Spnumber", groupvars=c("Time","Treatment")) Hbsp$Time = gsub("two weeks", "2", Hbsp$Time) Hbsp$Time = gsub("seven weeks", "7", Hbsp$Time) Hbsp$Time = gsub("twelve weeks","12", Hbsp$Time) Hbsp$Time = as.factor(Hbsp$Time) Hbsp$Time = ordered(Hbsp$Time,levels = c("2", "7", "12")) Hbsp$Treatment= as.factor(Hbsp$Treatment) Hbsp$Treatment <- ordered(Hbsp$Treatment,levels = c("WBB", "OAB", "OGB", "OBB")) summary(Hbsp) levels(Hbsp$Time) label_Hbsp <- c("", "", "", "", "a", "b", "a", "ab", "", "","", "") Hbsp_fig <- ggplot(Hbsp, aes(x=Time, y=Spnumber,color=Treatment, fill=Treatment)) + geom_bar(position=position_dodge(), stat="identity")+ geom_errorbar(aes(ymin=Spnumber-se, ymax=Spnumber+se), colour="black", width=.1, position=position_dodge(.9)) + geom_text(aes(y = Spnumber + 2 * se, label = label_Hbsp, group = Treatment), colour="red", position = position_dodge(0.8), size = 3, fontface = "bold") + ggtitle("C") + xlab("Sampling time (weeks)") + ylab(" species number per soil block") + scale_y_continuous(expand = c(0, 0)) + geom_text(aes(y = 12), label = "") + theme_bw() Hbsp_fig grid.newpage() pushViewport(viewport(layout = grid.layout(1, 3))) vplayout = function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) print(Hasp_fig, vp = vplayout(1, 1)) print(Hgsp_fig, vp = vplayout(1, 2)) print(Hbsp_fig, vp = vplayout(1, 3)) # PERMANOVA in Arable land (A) ------------------------------------------ library(vegan) HaT1WAOA = HaT1[c(1:5,6:10),6:19] HaT1WAOAgp = HaT1[c(1:5,6:10),c(1,3)] adonis(HaT1WAOA ~ Treatment, HaT1WAOAgp, permutations = 999, method="bray") HaT1WAOG = HaT1[c(1:5,11:15),6:19] HaT1WAOGgp= HaT1[c(1:5,11:15),c(1,3)] adonis(HaT1WAOG ~ Treatment, HaT1WAOGgp, permutations = 999, method="bray") HaT1WAOB = HaT1[c(1:5,16:20),6:19] HaT1WAOBgp= HaT1[c(1:5,16:20),c(1,3)] adonis(HaT1WAOB ~ Treatment, HaT1WAOBgp, permutations = 999, method="bray") HaT2WAOA = HaT2[c(1:5,6:10),6:19] HaT2WAOAgp= HaT2[c(1:5,6:10),c(1,3)] adonis(HaT2WAOA ~ Treatment, HaT2WAOAgp, permutations = 999, method="bray") HaT2WAOG = HaT2[c(1:5,11:15),6:19] HaT2WAOGgp= HaT2[c(1:5,11:15),c(1,3)] adonis(HaT2WAOG ~ Treatment, HaT2WAOGgp, permutations = 999, method="bray") HaT2WAOB = HaT2[c(1:5,16:20),6:19] HaT2WAOBgp= HaT2[c(1:5,16:20),c(1,3)] adonis(HaT2WAOB ~ Treatment, HaT2WAOBgp, permutations = 999, method="bray") HaT3WAOA = HaT3[c(1:5,6:10),6:19] HaT3WAOAgp= HaT3[c(1:5,6:10),c(1,3)] adonis(HaT3WAOA ~ Treatment, HaT3WAOAgp, permutations = 999, method="bray") HaT3WAOG = HaT3[c(1:5,11:15),6:19] HaT3WAOGgp= HaT3[c(1:5,11:15),c(1,3)] adonis(HaT3WAOG ~ Treatment, HaT3WAOGgp, permutations = 999, method="bray") HaT3WAOB = HaT3[c(1:5,16:20),6:19] HaT3WAOBgp= HaT3[c(1:5,16:20),c(1,3)] adonis(HaT3WAOB ~ Treatment, HaT3WAOBgp, permutations = 999, method="bray") # PERMANOVA in Grassland (G) ---------------------------------------------------------- HgT1WGOG = HgT1[c(1:5,6:10),6:19] HgT1WGOGgp= HgT1[c(1:5,6:10),c(1,3)] adonis(HgT1WGOG ~ Treatment, HgT1WGOGgp, permutations = 999, method="bray") HgT1WGOA = HgT1[c(1:5,11:15),6:19] HgT1WGOAgp= HgT1[c(1:5,11:15),c(1,3)] adonis(HgT1WGOA ~ Treatment, HgT1WGOAgp, permutations = 999, method="bray") HgT1WGOB = HgT1[c(1:5,16:20),6:19] HgT1WGOBgp= HgT1[c(1:5,16:20),c(1,3)] adonis(HgT1WGOB ~ Treatment, HgT1WGOBgp, permutations = 999, method="bray") HgT2WGOG = HgT2[c(1:5,6:10),6:19] HgT2WGOGgp= HgT2[c(1:5,6:10),c(1,3)] adonis(HgT2WGOG ~ Treatment, HgT2WGOGgp, permutations = 999, method="bray") HgT2WGOA = HgT2[c(1:5,11:15),6:19] HgT2WGOAgp= HgT2[c(1:5,11:15),c(1,3)] adonis(HgT2WGOA ~ Treatment, HgT2WGOAgp, permutations = 999, method="bray") HgT2WGOB = HgT2[c(1:5,16:20),6:19] HgT2WGOBgp= HgT2[c(1:5,16:20),c(1,3)] adonis(HgT2WGOB ~ Treatment, HgT2WGOBgp, permutations = 999, method="bray") HgT3WGOG = HgT3[c(1:5,6:10),6:19] HgT3WGOGgp= HgT3[c(1:5,6:10),c(1,3)] adonis(HgT3WGOG ~ Treatment, HgT3WGOGgp, permutations = 999, method="bray") HgT3WGOA = HgT3[c(1:5,11:15),6:19] HgT3WGOAgp= HgT3[c(1:5,11:15),c(1,3)] adonis(HgT3WGOA ~ Treatment, HgT3WGOAgp, permutations = 999, method="bray") HgT3WGOB = HgT3[c(1:5,16:20),6:19] HgT3WGOBgp= HgT3[c(1:5,16:20),c(1,3)] adonis(HgT3WGOB ~ Treatment, HgT3WGOBgp, permutations = 999, method="bray") # PERMANOVA in BSAL ---------------------------------------------------------- HbT1WBOB = HbT1[c(1:5,6,8:10),6:19] HbT1WBOBgp= HbT1[c(1:5,6,8:10),c(1,3)] adonis(HbT1WBOB ~ Treatment, HbT1WBOBgp, permutations = 999, method="bray") HbT1WBOA = HbT1[c(1:5,11:15),6:19] HbT1WBOAgp= HbT1[c(1:5,11:15),c(1,3)] adonis(HbT1WBOA ~ Treatment, HbT1WBOAgp, permutations = 999, method="bray") HbT1WBOG = HbT1[c(1:5,16:20),6:19] HbT1WBOGgp= HbT1[c(1:5,16:20),c(1,3)] adonis(HbT1WBOG ~ Treatment, HbT1WBOGgp, permutations = 999, method="bray") HbT2WBOB = HbT2[c(1:5,6:8,10),6:19] HbT2WBOBgp= HbT2[c(1:5,6:8,10),c(1,3)] adonis(HbT2WBOB ~ Treatment, HbT2WBOBgp, permutations = 999, method="bray") HbT2WBOA = HbT2[c(1:5,11:15),6:19] HbT2WBOAgp= HbT2[c(1:5,11:15),c(1,3)] adonis(HbT2WBOA ~ Treatment, HbT2WBOAgp, permutations = 999, method="bray") HbT2WBOG = HbT2[c(1:5,16:20),6:19] HbT2WBOGgp= HbT2[c(1:5,16:20),c(1,3)] adonis(HbT2WBOG ~ Treatment, HbT2WBOGgp, permutations = 999, method="bray") HbT3WBOB = HbT3[c(1:5,6:10),6:19] HbT3WBOBgp= HbT3[c(1:5,6:10),c(1,3)] adonis(HbT3WBOB ~ Treatment, HbT3WBOBgp, permutations = 999, method="bray") HbT3WBOA = HbT3[c(1:5,11:15),6:19] HbT3WBOAgp= HbT3[c(1:5,11:15),c(1,3)] adonis(HbT3WBOA ~ Treatment, HbT3WBOAgp, permutations = 999, method="bray") HbT3WBOG = HbT3[c(1:5,16:20),6:19] HbT3WBOGgp= HbT3[c(1:5,16:20),c(1,3)] adonis(HbT3WBOG ~ Treatment, HbT3WBOGgp, permutations = 999, method="bray") #PCoA analyses library(vegan) library(ggplot2) library(plyr) #input species table and group information recoll <- read.csv("Coll_recov_all180_nz.csv", row.names = 1, header = TRUE) recoll2 <- read.csv("Coll_recov_all180_nz.csv", header = TRUE) #in the Arabland (A), 2 weeks recoll.a2 <- subset(recoll, Destination=="A" & Time == "two weeks") colSums(recoll.a2[,5:18]) rowSums(recoll.a2[,5:18]) #Keep species with more than one individuals in all samples sp_tb.a2 <- recoll.a2[,c(5:7,9,12:13,15:17)] recoll2.a2 <-subset(recoll2, Destination=="A" & Time == "two weeks") group.a2 <- recoll2.a2[,1:4] #PCoA Ordination of A, 2 weeks distance.a2 <- vegdist(sp_tb.a2, method = 'bray') pcoa.a2 <- cmdscale(distance.a2, k = (nrow(sp_tb.a2) - 1), eig = T) summary(pcoa.a2) #ggplot2 PCoA for A, 2 weeks pcoa.a2_eig <- (pcoa.a2$eig)[1:2] / sum(pcoa.a2$eig) sample_site.a2 <- data.frame({pcoa.a2$point})[1:2] sample_site.a2$name <- rownames(sample_site.a2) names(sample_site.a2)[1:2] <- c('PCoA1', 'PCoA2') sample_site.a2 <- merge(sample_site.a2, group.a2, by = 'name', all.x = T) sample_site.a2$Treatment <- factor(sample_site.a2$Treatment, levels = c('WAA', 'OAA', 'OGA', 'OBA')) group.a2_border <- ddply(sample_site.a2, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.a2_plot <- ggplot(sample_site.a2, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.a2_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(15,0,12,7)) + scale_fill_manual(values = c('#0000FF2E', '#008B002E','#9400D32E','#FFA5002E')) + guides(shape = guide_legend(order = 1),fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.a2_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.a2_eig[2], 2), '%')) pcoa.a2_plot #Arabland (A), 7 weeks recoll.a7 <- subset(recoll, Destination=="A" & Time == "seven weeks") colSums(recoll.a7[,5:18]) rowSums(recoll.a7[,5:18]) sp_tb.a7 <- recoll.a7[,c(5:9,12,14:18)] recoll2.a7 <-subset(recoll2, Destination=="A" & Time == "seven weeks") group.a7 <- recoll2.a7[,1:4] distance.a7 <- vegdist(sp_tb.a7, method = 'bray') pcoa.a7 <- cmdscale(distance.a7, k = (nrow(sp_tb.a7) - 1), eig = T) summary(pcoa.a7) pcoa.a7_eig <- (pcoa.a7$eig)[1:2] / sum(pcoa.a7$eig) sample_site.a7 <- data.frame({pcoa.a7$point})[1:2] sample_site.a7$name <- rownames(sample_site.a7) names(sample_site.a7)[1:2] <- c('PCoA1', 'PCoA2') sample_site.a7 <- merge(sample_site.a7, group.a7, by = 'name', all.x = T) sample_site.a7$Treatment <- factor(sample_site.a7$Treatment, levels = c('WAA', 'OAA', 'OGA', 'OBA')) group.a7_border <- ddply(sample_site.a7, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.a7_plot <- ggplot(sample_site.a7, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.a7_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(15,0,12,7)) + scale_fill_manual(values = c('#0000FF2E', '#008B002E','#9400D32E','#FFA5002E')) + guides(shape = guide_legend(order = 1), fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.a7_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.a7_eig[2], 2), '%')) pcoa.a7_plot #Arabland (A), 12 weeks recoll.a12 <- subset(recoll, Destination=="A" & Time == "twelve weeks") colSums(recoll.a12[,5:18]) rowSums(recoll.a12[,5:18]) sp_tb.a12 <- recoll.a12[,c(5:7,9:10,12:18)] recoll2.a12 <-subset(recoll2, Destination=="A" & Time == "twelve weeks") group.a12 <- recoll2.a12[,1:4] distance.a12 <- vegdist(sp_tb.a12, method = 'bray') pcoa.a12 <- cmdscale(distance.a12, k = (nrow(sp_tb.a12) - 1), eig = T) summary(pcoa.a12) pcoa.a12_eig <- (pcoa.a12$eig)[1:2] / sum(pcoa.a12$eig) sample_site.a12 <- data.frame({pcoa.a12$point})[1:2] sample_site.a12$name <- rownames(sample_site.a12) names(sample_site.a12)[1:2] <- c('PCoA1', 'PCoA2') sample_site.a12 <- merge(sample_site.a12, group.a12, by = 'name', all.x = T) sample_site.a12$Treatment <- factor(sample_site.a12$Treatment, levels = c('WAA', 'OAA', 'OGA', 'OBA')) group.a12_border <- ddply(sample_site.a12, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.a12_plot <- ggplot(sample_site.a12, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.a12_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(15,0,12,7)) + scale_fill_manual(values = c('#0000FF2E', '#008B002E','#9400D32E','#FFA5002E')) + guides(shape = guide_legend(order = 1), fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.a12_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.a12_eig[2], 2), '%')) pcoa.a12_plot #Grassland (G), 2 weeks recoll.g2 <- subset(recoll, Destination=="G" & Time == "two weeks") colSums(recoll.g2[,5:18]) rowSums(recoll.g2[,5:18]) sp_tb.g2 <- recoll.g2[,c(5:7,9,11:12,15:17)] recoll2.g2 <-subset(recoll2, Destination=="G" & Time == "two weeks") group.g2 <- recoll2.g2[,1:4] distance.g2 <- vegdist(sp_tb.g2, method = 'bray') pcoa.g2 <- cmdscale(distance.g2, k = (nrow(sp_tb.g2) - 1), eig = T) summary(pcoa.g2) pcoa.g2_eig <- (pcoa.g2$eig)[1:2] / sum(pcoa.g2$eig) sample_site.g2 <- data.frame({pcoa.g2$point})[1:2] sample_site.g2$name <- rownames(sample_site.g2) names(sample_site.g2)[1:2] <- c('PCoA1', 'PCoA2') sample_site.g2 <- merge(sample_site.g2, group.g2, by = 'name', all.x = T) sample_site.g2$Treatment <- factor(sample_site.g2$Treatment, levels = c('WGG', 'OGG', 'OAG', 'OBG')) group.g2_border <- ddply(sample_site.g2, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.g2_plot <- ggplot(sample_site.g2, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.g2_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(17,2,5,9)) + scale_fill_manual(values = c('#0000FF2E','#9400D32E', '#008B002E','#FFA5002E')) + guides(shape = guide_legend(order = 1), fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.g2_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.g2_eig[2], 2), '%')) pcoa.g2_plot #Grassland (G), 7 weeks recoll.g7 <- subset(recoll, Destination=="G" & Time == "seven weeks") colSums(recoll.g7[,5:18]) rowSums(recoll.g7[,5:18]) sp_tb.g7 <- recoll.g7[,c(5,7,12:16,18)] recoll2.g7 <-subset(recoll2, Destination=="G" & Time == "seven weeks") group.g7 <- recoll2.g7[,1:4] distance.g7 <- vegdist(sp_tb.g7, method = 'bray') pcoa.g7 <- cmdscale(distance.g7, k = (nrow(sp_tb.g7) - 1), eig = T) summary(pcoa.g7) pcoa.g7_eig <- (pcoa.g7$eig)[1:2] / sum(pcoa.g7$eig) sample_site.g7 <- data.frame({pcoa.g7$point})[1:2] sample_site.g7$name <- rownames(sample_site.g7) names(sample_site.g7)[1:2] <- c('PCoA1', 'PCoA2') sample_site.g7 <- merge(sample_site.g7, group.g7, by = 'name', all.x = T) sample_site.g7$Treatment <- factor(sample_site.g7$Treatment, levels = c('WGG', 'OGG', 'OAG', 'OBG')) group.g7_border <- ddply(sample_site.g7, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.g7_plot <- ggplot(sample_site.g7, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.g7_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(17,2,5,9)) + scale_fill_manual(values = c('#0000FF2E','#9400D32E', '#008B002E','#FFA5002E')) + guides(shape = guide_legend(order = 1), fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.g7_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.g7_eig[2], 2), '%')) pcoa.g7_plot #Grassland (G), 12 weeks recoll.g12 <- subset(recoll, Destination=="G" & Time == "twelve weeks") colSums(recoll.g12[,5:18]) rowSums(recoll.g12[,5:18]) sp_tb.g12 <- recoll.g12[,c(5,7,9,12:14,16:17)] recoll2.g12 <-subset(recoll2, Destination=="G" & Time == "twelve weeks") group.g12 <- recoll2.g12[,1:4] distance.g12 <- vegdist(sp_tb.g12, method = 'bray') pcoa.g12 <- cmdscale(distance.g12, k = (nrow(sp_tb.g12) - 1), eig = T) summary(pcoa.g12) pcoa.g12_eig <- (pcoa.g12$eig)[1:2] / sum(pcoa.g12$eig) sample_site.g12 <- data.frame({pcoa.g12$point})[1:2] sample_site.g12$name <- rownames(sample_site.g12) names(sample_site.g12)[1:2] <- c('PCoA1', 'PCoA2') sample_site.g12 <- merge(sample_site.g12, group.g12, by = 'name', all.x = T) sample_site.g12$Treatment <- factor(sample_site.g12$Treatment, levels = c('WGG', 'OGG', 'OAG', 'OBG')) group.g12_border <- ddply(sample_site.g12, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.g12_plot <- ggplot(sample_site.g12, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.g12_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(17,2,5,9)) + scale_fill_manual(values = c('#0000FF2E','#9400D32E', '#008B002E','#FFA5002E')) + guides(shape = guide_legend(order = 1), fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.g12_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.g12_eig[2], 2), '%')) pcoa.g12_plot #BSAL (B), 2 weeks recoll.b2 <- subset(recoll, Destination=="B" & Time == "two weeks") colSums(recoll.b2[,5:18]) rowSums(recoll.b2[,5:18]) sp_tb.b2 <- recoll.b2[-7,c(5:7,10,12,15:17)] recoll2.b2 <-subset(recoll2, Destination=="B" & Time == "two weeks") group.b2 <- recoll2.b2[,1:4] distance.b2 <- vegdist(sp_tb.b2, method = 'bray') pcoa.b2 <- cmdscale(distance.b2, k = (nrow(sp_tb.b2) - 1), eig = T) summary(pcoa.b2) pcoa.b2_eig <- (pcoa.b2$eig)[1:2] / sum(pcoa.b2$eig) sample_site.b2 <- data.frame({pcoa.b2$point})[1:2] sample_site.b2$name <- rownames(sample_site.b2) names(sample_site.b2)[1:2] <- c('PCoA1', 'PCoA2') sample_site.b2 <- merge(sample_site.b2, group.b2, by = 'name', all.x = T) sample_site.b2$Treatment <- factor(sample_site.b2$Treatment, levels = c('WBB', 'OBB', 'OAB', 'OGB')) group.b2_border <- ddply(sample_site.b2, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.b2_plot <- ggplot(sample_site.b2, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.b2_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(19, 1, 10, 13)) + scale_fill_manual(values = c('#0000FF2E','#FFA5002E', '#008B002E','#9400D32E')) + guides(shape = guide_legend(order = 1), fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.b2_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.b2_eig[2], 2), '%')) pcoa.b2_plot #BSAL (B), 7 weeks recoll.b7 <- subset(recoll, Destination=="B" & Time == "seven weeks") colSums(recoll.b7[,5:18]) rowSums(recoll.b7[,5:18]) sp_tb.b7 <- recoll.b7[c(-8,-9),c(5:7,12,15:16)] recoll2.b7 <-subset(recoll2, Destination=="B" & Time == "seven weeks") group.b7 <- recoll2.b7[,1:4] distance.b7 <- vegdist(sp_tb.b7, method = 'bray') pcoa.b7 <- cmdscale(distance.b7, k = (nrow(sp_tb.b7) - 1), eig = T) summary(pcoa.b7) pcoa.b7_eig <- (pcoa.b7$eig)[1:2] / sum(pcoa.b7$eig) sample_site.b7 <- data.frame({pcoa.b7$point})[1:2] sample_site.b7$name <- rownames(sample_site.b7) names(sample_site.b7)[1:2] <- c('PCoA1', 'PCoA2') sample_site.b7 <- merge(sample_site.b7, group.b7, by = 'name', all.x = T) sample_site.b7$Treatment <- factor(sample_site.b7$Treatment, levels = c('WBB', 'OBB', 'OAB', 'OGB')) group.b7_border <- ddply(sample_site.b7, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.b7_plot <- ggplot(sample_site.b7, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.b7_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(19, 1, 10, 13)) + scale_fill_manual(values = c('#0000FF2E','#FFA5002E', '#008B002E','#9400D32E')) + guides(shape = guide_legend(order = 1), fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.b7_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.b7_eig[2], 2), '%')) pcoa.b7_plot #BSAL (B), 12 weeks recoll.b12 <- subset(recoll, Destination=="B" & Time == "twelve weeks") colSums(recoll.b12[,5:18]) rowSums(recoll.b12[,5:18]) sp_tb.b12 <- recoll.b12[,c(5,7,12,16,18)] recoll2.b12 <-subset(recoll2, Destination=="B" & Time == "twelve weeks") group.b12 <- recoll2.b12[,1:4] distance.b12 <- vegdist(sp_tb.b12, method = 'bray') pcoa.b12 <- cmdscale(distance.b12, k = (nrow(sp_tb.b12) - 1), eig = T) summary(pcoa.b12) pcoa.b12_eig <- (pcoa.b12$eig)[1:2] / sum(pcoa.b12$eig) sample_site.b12 <- data.frame({pcoa.b12$point})[1:2] sample_site.b12$name <- rownames(sample_site.b12) names(sample_site.b12)[1:2] <- c('PCoA1', 'PCoA2') sample_site.b12 <- merge(sample_site.b12, group.b12, by = 'name', all.x = T) sample_site.b12$Treatment <- factor(sample_site.b12$Treatment, levels = c('WBB', 'OBB', 'OAB', 'OGB')) group.b12_border <- ddply(sample_site.b12, 'Treatment', function(df) df[chull(df[[2]], df[[3]]), ]) pcoa.b12_plot <- ggplot(sample_site.b12, aes(PCoA1, PCoA2, group = Treatment)) + theme(panel.grid = element_line(color = 'gray', linetype = 2, size = 0.1), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) + geom_vline(xintercept = 0, color = 'gray', size = 0.4) + geom_hline(yintercept = 0, color = 'gray', size = 0.4) + geom_polygon(data = group.b12_border, aes(fill = Treatment)) + geom_point(aes(shape = Treatment), size = 2.5, alpha = 0.8) + scale_shape_manual(values = c(19, 1, 10, 13)) + scale_fill_manual(values = c('#0000FF2E','#FFA5002E', '#008B002E','#9400D32E')) + guides(shape = guide_legend(order = 1), fill = guide_legend(order = 2)) + labs(x = paste('PCoA axis1: ', round(100 * pcoa.b12_eig[1], 2), '%'), y = paste('PCoA axis2: ', round(100 * pcoa.b12_eig[2], 2), '%')) pcoa.b12_plot library(grid) grid.newpage() pushViewport(viewport(layout = grid.layout(3, 3))) vplayout = function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) print(pcoa.a2_plot, vp = vplayout(1, 1)) print(pcoa.a7_plot, vp = vplayout(1, 2)) print(pcoa.a12_plot, vp = vplayout(1, 3)) print(pcoa.g2_plot, vp = vplayout(2, 1)) print(pcoa.g7_plot, vp = vplayout(2, 2)) print(pcoa.g12_plot, vp = vplayout(2, 3)) print(pcoa.b2_plot, vp = vplayout(3, 1)) print(pcoa.b7_plot, vp = vplayout(3, 2)) print(pcoa.b12_plot, vp = vplayout(3, 3))