#### Description: This script replicates Fig. 1 and Fig. S6 in the manuscript, "Multiple plant traits shape the genetic basis of herbivore community assembly." source('~/Documents/Genotype_Willow_Community/datasets_&_Rscripts/arthropods_&_traits/herb_trait_manage.r') ## load libraries for plotting library(ggplot2) library(gridExtra) library(GGally) library(extrafont) library(extrafontdb) library(Rttf2pt1) ## melt herbivore community data library(reshape2) melt.com <- melt(herb.com, id.vars=c("Gender","Genotype","Plant.Position")) ## create basic attributes for plots points <- geom_point(size = 5, shape = 21, fill = "white") errors <- geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = 0.2) theme <- theme_bw() + theme(text = element_text(family = "Arial", size = 18), legend.key = element_blank(), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 0.5), axis.line=element_line(colour="black"), panel.grid=element_blank(), panel.border=element_blank(), legend.title=element_blank() ) ## Note that both the herbivore abundance and richness plots contain the mean and SE of each responses for each willow genotype. Genotypes were ordered based on mean herbivore abundance. Also, Genotype '*' was changed to 'C' for plotting aesethetics. library(plyr) ## Total herbivore abundance plot (Fig. 1B) herbabund.geno.sum <- ddply(ddply(subset(melt.com, variable == "Herbivore_abundance"), .(Genotype,variable), summarise, mean=mean(value), se = sd(value)/sqrt(length(value))), .(mean)) # calculates mean and se herbabund.geno.order <- reorder(herbabund.geno.sum$Genotype, herbabund.geno.sum$mean) levels(herbabund.geno.order)[15] <- "C" # '*' changed to "C" for aesthetic purposes (herbabund.plot <- ggplot(herbabund.geno.sum, aes(x = herbabund.geno.order, y = mean)) + errors + xlab(expression(paste(italic("S. hookeriana")," genotype"))) + ylab("No. of individuals") + points + theme + coord_cartesian(ylim=c(19,101)) + scale_y_continuous(breaks = seq(20,110,20))) ## Herbivore richness plot (Fig. 1A) herbrich.geno.sum <- ddply(ddply(subset(melt.com, variable == "herbivore_rich"), .(Genotype,variable), summarise, mean=mean(value), se = sd(value)/sqrt(length(value))), .(mean)) herbrich.geno.order <- factor(herbrich.geno.sum$Genotype, levels = levels(herbabund.geno.order)) (herbrich.plot <- ggplot(herbrich.geno.sum, aes(x = herbrich.geno.order, y = mean)) + errors + xlab("") + ylab("No. of species") + points + theme + theme(axis.text.x = element_blank())+ coord_cartesian(ylim=c(8,20)) + scale_y_continuous(breaks = seq(8,20,2)) ) ## resize plots for a multipanelled figure. Saved plot as a pdf file and exported as a 5x5 figure. gp1<- ggplot_gtable(ggplot_build(herbrich.plot)) gp2<- ggplot_gtable(ggplot_build(herbabund.plot)) maxWidth = unit.pmax(gp1$widths[2:3], gp2$widths[2:3]) gp1$widths[2:3] <- maxWidth gp2$widths[2:3] <- maxWidth grid.arrange(gp1, gp2, ncol = 1) #### Create herbivore community composition plot (Fig. 1c). Note that I created this figure separately and joined it with the previous figure in a different program to make sure all of the axes were aligned. ## perform RDA analyses and extract centroid scores for plotting. library(vegan) comp.rda.geno.only <- rda(decostand(herb.com[ ,herbivore.names], method = "normalize") ~ Genotype, data = herb.com) summary(comp.rda.geno.only) centroids.herb.geno <- scores(comp.rda.geno.only, choices=c(1,2), display="cn") row.names(centroids.herb.geno) <- c("C","A","B", LETTERS[4:26]) # switches genotype "C" for "*" for aesthetic purposes only. then preserves the correct order. ## generate plot. Coded out the manner in which I saved the plot # pdf(file="herb.com.plot.pdf",width=5.5, height=5.5) plot.new() par(family="Arial") plot.window( xlim=c(-0.55,0.58), ylim = c(-0.69,0.5), asp = 1) abline(h=0, lty="dotted") abline(v=0, lty="dotted") axis(side=1, tck = -0.015, padj=-0.5) axis(side=2, tck= -0.015, las=1) title(xlab="RDA 1 (7.4%)", ylab="RDA 2 (5.1%)") box() # fill in the plot ordiellipse(comp.rda.geno.only, groups = herb.com$Genotype, kind = "se", draw = "lines", col="grey") ordiellipse(comp.rda.geno.only, groups = herb.com$Genotype, kind = "se", show.groups="E", col = "grey", draw = "polygon") ordiellipse(comp.rda.geno.only, groups = herb.com$Genotype, kind = "se", show.groups="*", col = "grey", draw = "polygon") # note the C and I switch for aesthetic reasons ordiellipse(comp.rda.geno.only, groups = herb.com$Genotype, kind = "se", show.groups="L", col = "grey", draw = "polygon") ordiellipse(comp.rda.geno.only, groups = herb.com$Genotype, kind = "se", show.groups="R", col = "grey", draw = "polygon") ordiellipse(comp.rda.geno.only, groups = herb.com$Genotype, kind = "se", show.groups="Z", col = "grey", draw = "polygon") text(x = centroids.herb.geno[ ,1], centroids.herb.geno[ ,2], labels = row.names(centroids.herb.geno)) #### Generate herbivore community composition plot including the focal herbivore species (Fig. S6) ## give genus of focal herbivores selecting.sp <- ordiplot(comp.rda.geno.only, display="species") row.names(selecting.sp[[1]])[c(1,2,8,9,17)] <- c("Iteomyia","Aculus","Caloptilia", "Tachyerges","Cicadellid nymph sp. 1") ## create blank plot # pdf(file="herb.com.plot.withSpecies.pdf",width=5.5, height=5.5) # code for generating pdf plot.new() plot.window( xlim=c(-0.55,0.58), ylim = c(-0.69,0.5), asp = 1) abline(h=0, lty="dotted") abline(v=0, lty="dotted") axis(side=1) axis(side=2) title(xlab="RDA 1 (7.4%)", ylab="RDA 2 (5.1%)") box() ## fill in the plot text(x = centroids.herb.geno[ ,1], centroids.herb.geno[ ,2], labels = row.names(centroids.herb.geno), col = "grey") points(comp.rda.geno.only, display="sp") identify(selecting.sp, 'species') # press 'esc' button on computer after manually selecting species in the plot. # dev.off() # turn off pdf device.