#### Description: This script runs one-way ANOVA, GLMs, and RDA models of herbivore community responses to willow genotype. These results are reported in the first 2 paragraphs of the Results section of the manuscript entitled, "Multiple plant traits shape the genetic basis of herbivore community assembly." ## source in required scripts and load up required libraries. source('~/Documents/Genotype_Willow_Community/datasets_&_Rscripts/arthropods_&_traits/herb_trait_manage.r') library(car) # for 'logit' transformation library(vegan) # for community composition analyses library(MASS) # for negative binomial GLM ## function for compiling data table anova_table <- function(name, y, model){ G = herb.com$Genotype min_max = data.frame(min = min(aggregate(y ~ G, FUN = mean)[,2]), max = max(aggregate(y ~ G, FUN = mean,)[,2])) Test = round(get(paste("model"))[1,4],2) P = round(get(paste("model"))[1,5],3) c(name, round(min_max[,1],1), round(min_max[,2],1), Test, P) } #### Run ANOVAs for herbivore abundance, richness, rarefied richness, and evenness. ## log(herbivore abundance) abund.anova <- anova(aov(log(Herbivore_abundance) ~ Genotype, herb.com)) abund.anova2 <- anova_table(name = "log(Herbivore abundance)", y = herb.com$Herbivore_abundance, model = abund.anova) ## log(herbivore richness) rich.anova <- anova(aov(log(herbivore_rich) ~ Genotype, herb.com)) rich.anova2 <- anova_table(name = "log(Herbivore richness)", y = herb.com$herbivore_rich, model = rich.anova) ## herbivore rarefied richness rrich.anova <- anova(aov(Herbivores_rarerich ~ Genotype, herb.com)) rrich.anova2 <- anova_table(name = "Herbivore rarefied richness", y = herb.com$Herbivores_rarerich, model = rrich.anova) ## logit(herbivore evenness) even.anova <- anova(aov(logit(Herbivores_evenness) ~ Genotype, herb.com)) even.anova2 <- anova_table(name = "logit(Herbivore evenness)", y = herb.com$Herbivores_evenness, model = even.anova) #### Run negative binomial GLMs of herbivore feeding guilds and species responses to willow genotype. ## leaf chewer damage (percent leaf area remvoed) ld.anova <- anova(aov(log(avg.LeafDam.score) ~ Genotype, herb.com)) ld.anova2 <- anova_table(name = "log(avg. leaf dam score)", y = herb.com$avg.LeafDam.score, model = ld.anova) ## Gallers gall0 <- glm.nb(gallers_abund ~ 1, herb.com) gall1 <- update(gall0, .~. + Genotype) gall.LR <- anova(gall0,gall1) gall.LR <- gall.LR[-1,-c(1:3)] # remove initial gall.LR <- anova_table(name = "Gall abund glm.nb", y = herb.com$gallers_abund, model = gall.LR) ## Leaf miners leafmine0 <- glm.nb(leaf_miner_abund ~ 1, herb.com) leafmine1 <- update(leafmine0, .~. + Genotype) leafmine.LR <- anova(leafmine0,leafmine1) leafmine.LR <- leafmine.LR[-1,-c(1:3)] # remove initial leafmine.LR <- anova_table(name = "Leaf miner abund glm.nb", y = herb.com$leaf_miner_abund, model = leafmine.LR) ## Xylem suckers xylemsuck0 <- glm.nb(xylem_suckers_abund ~ 1, herb.com) xylemsuck1 <- update(xylemsuck0, .~. + Genotype) xylemsuck.LR <- anova(xylemsuck0,xylemsuck1) xylemsuck.LR <- xylemsuck.LR[-1,-c(1:3)] # remove initial xylemsuck.LR <- anova_table(name = "xylem feeder abund glm.nb", y = herb.com$xylem_suckers_abund, model = xylemsuck.LR) ## Phloem suckers. phloemsuck0 <- glm.nb(phloem_suckers_abund ~ 1, herb.com) phloemsuck1 <- update(phloemsuck0, .~. + Genotype) phloemsuck.LR <- anova(phloemsuck0,phloemsuck1) phloemsuck.LR <- phloemsuck.LR[-1,-c(1:3)] # remove initial phloemsuck.LR <- anova_table(name = "Phloem feeder abund glm.nb", y = herb.com$phloem_suckers_abund, model = phloemsuck.LR) ## Leaf chewer abundance. leafchew0 <- glm.nb(leaf_chewers_total_abund ~ 1, herb.com) leafchew1 <- update(leafchew0, .~. + Genotype) leafchew.LR <- anova(leafchew0,leafchew1) leafchew.LR <- leafchew.LR[-1,-c(1:3)] # remove initial leafchew.LR <- anova_table(name = "Leaf chewer abund glm.nb", y = herb.com$leaf_chewers_total_abund, model = leafchew.LR) #### Redundancy analysis (RDA) of herbivore community composition among willow genotypes. Used chord transformation of community data to focus on differences in relative composition of herbivore species among genotypes. ## First, check whether there is overdispersion. If there isn't, then there is more support for a mean difference in centroids of community composition comp.betadisper <- betadisper(vegdist(decostand(herb.com[ ,herbivore.names], method = "normalize"), method = "euclidean"), herb.com$Genotype, bias.adjust=TRUE) anova(comp.betadisper) # no evidence of overdispersion ## Redundancy Analysis (RDA) comp.rda <- rda(decostand(herb.com[ ,herbivore.names], method = "normalize") ~ Genotype, data = herb.com) summary(comp.rda) comp.rda.anova <- anova(comp.rda, step=1000) comp.rda.anova <- c("Community composition RDA", "","",round(comp.rda.anova[1,3],2), round(comp.rda.anova[1,5],3)) plot(comp.rda, display = "sp") # identify species driving variation among Genotypes # mLG = Aculups tetanothrix # LTF = Caloptilia sp. # vLG = Iteomyia salicisverruca # Cicadellidae nymph 1 = Cicadellidae sp. 1 # Curculionid_1 = Tachyerges salicis #### Negative binomial GLMs of Individual Herbivore responses to willow genotype. These analyses were restricted to the herbivores that drove the variation in community response. ## Cicadellidae nymph. cicnymph0 <- glm.nb(Cicadellidae_nymph_1 ~ 1, herb.com) cicnymph1 <- update(cicnymph0, .~. + Genotype) cicnymph.LR <- anova(cicnymph0,cicnymph1) cicnymph.LR <- cicnymph.LR[-1,-c(1:3)] # remove initial cicnymph.LR <- anova_table(name = "Cicadellid nymph sp. 1 abund glm.nb", y = herb.com$Cicadellidae_nymph_1, model = cicnymph.LR) ## Curculionid_1 (Tachyerges salicis). curc10 <- glm.nb(Curculionid_1 ~ 1, herb.com) curc11 <- update(curc10, .~. + Genotype) curc1.LR <- anova(curc10,curc11) curc1.LR <- curc1.LR[-1,-c(1:3)] # remove initial curc1.LR <- anova_table(name = "Tachyerges salicis abund glm.nb", y = herb.com$Curculionid_1, model = curc1.LR) ## LTF = Caloptilia sp. LTF0 <- glm.nb(LTF ~ 1, herb.com) LTF1 <- update(LTF0, .~. + Genotype) LTF.LR <- anova(LTF0,LTF1) LTF.LR <- LTF.LR[-1,-c(1:3)] # remove initial LTF.LR <- anova_table(name = "Caloptilia sp. abund glm.nb", y = herb.com$LTF, model = LTF.LR) ## vLG = Iteomyia salicisverruca. vLG0 <- glm.nb(vLG ~ 1, herb.com) vLG1 <- update(vLG0, .~. + Genotype) vLG.LR <- anova(vLG0,vLG1) vLG.LR <- vLG.LR[-1,-c(1:3)] # remove initial vLG.LR <- anova_table(name = "Iteomyia salicisverruca abund glm.nb", y = herb.com$vLG, model = vLG.LR) Q <- which(herb.com$Genotype == "Q") # used to examine whether removing Genotype Q, which never had vLG, drove the effect of genotype (it wasn't). ## mLG = Aculus tetanothrix (mite leaf galler). mLG0 <- glm.nb(mLG ~ 1, herb.com) mLG1 <- update(mLG0, .~. + Genotype) mLG.LR <- anova(mLG0,mLG1) mLG.LR <- mLG.LR[-1,-c(1:3)] # remove initial mLG.LR <- anova_table(name = "Aculups tetanothrix abund glm.nb", y = herb.com$mLG, model = mLG.LR) ## Compile results from all herbivore responses to willow genotype into a dataframe anovas.data <- rbind.data.frame(rich.anova2, abund.anova2, rrich.anova2, even.anova2, comp.rda.anova, curc1.LR, LTF.LR, vLG.LR, mLG.LR, cicnymph.LR, ld.anova2, leafchew.LR, xylemsuck.LR, phloemsuck.LR, gall.LR, leafmine.LR) colnames(anovas.data) = c("Response","Min","Max","Test","P")