knitr::opts_chunk$set(echo = TRUE, out.width = '100%')
setwd("~/PhD")

Description of research

This study seeks to examine the degree of universality in plant defenses by comparing mass gain from a generalist herbivore in laboratory assays with naturally occurring arthropod counts on native plants in a common garden.

The metadata used for our analyses was collected annually in April-May for three consecutive years from 2021-2023 in a common garden at Concordia University Irvine. The common garden consisted of 14 species of woody shrubs native to the Coastal Sage Scrub ecosystem, arranged in 20 plots, each with 1 individual plant per species.

Data collected, abbreviations used and variable type assigned for the metadata can be found in the table below.

Data Description Column Name Variable Type
Plant # Individual plant Plant Factor
Plant Species Plant Species Species Factor
Plot Individual Plot Plot Factor
Year Year data collected Year Factor
Average Spodoptera mass (mg) Average mass of surviving Spodoptera Avg_spod_mass_mg Numeric
Total arthropod count Total number of arthropods Abundance_Total Numeric
Herbivore count Total number of herbivores Abundance_Herb Numeric

Methods

Spodoptera mass was collected after a 10-day laboratory feeding assay. Briefly, Spodoptera exigua, a generalist herbivore, were hatched from eggs and immediately fed leaves from each individual plant. At the end of 10 days, the number of surviving larvae and final mass of each larvae were recorded. See You, An, and Li (2020) for details. Plant and arthropod data were collected from each surviving plant from 2021-2023. Plant estimated biomass was collected using a branch-biomass estimation method as previously described in Nell and Mooney (2019). Arthropods were collected from plants using a vacuum method as previously described in Pratt et al. (2017), sorted from plant material and identified to taxonomic order.

1- Load packages and data tables

1.1 Load packages and data tables

library(dplyr)
library(lme4)
library(RColorBrewer)
library(devtools)
library(sjPlot)
library(tidyr)
library(broom.mixed)
library(knitr)
library(broom)
library(jtools)
library(kableExtra)
library(huxtable)
library(gtsummary)
library(car)
library(lmerTest)
library(ggplot2)
library(sjmisc)
library(sjstats)
library(sjlabelled)
library(glmmTMB)
library(DHARMa)
library(viridisLite)
library(data.table)
library(effects)
library(emmeans)
library(kableExtra)
library(MuMIn)
library(performance)
library(gridExtra)
library(cowplot)

2- Visualize raw data

Plot histograms of Herbivore abundance and Avg. Spodoptera mass to visualize distribution.

TAKEAWAY- Strongly skewed distributions. Log transformation will be applied to normalize data.

3- Community-level analyses

These sets of analyses will examine the relationship between Herbivore density on individual plants and Spodoptera growth rate across years.

Analyses will first be conducted at the level of the ecological community to examine broad-scale relationships between Spodoptera growth rate and herbivore abundance irrespective of Species-level differences. As such, Species, along with Plot, and Individual Plant, will included as random factors in the analysis. Year will be added as a fixed effect to account for year-year variation.

The generalized linear mixed model (lmer) that we will use to fit the data is:

Herbivore density ~ Spodoptera growth rate * Year + Spodoptera growth rate + Year + (1|Species) + (1|Plot) + (1|Plant)

3.1- ANOVA of Herbivore count ~ Spodoptera growth rate

We will examine the relationship between Spodoptera growth rate and Herbivore density across years while controlling for the effects of species, plant and plot on this relationship.

ANOVA results of log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Marginal_R2
log(Avg_spod_mass_mg) 0.5173969 1 0.4719542 0.07650765
Year 34.7109562 2 0.0000000 0.05354543
log(Avg_spod_mass_mg):Year 8.7604797 2 0.0125224 0.16071688

RESULTS

  • NO significance in relationship between Spodoptera growth rate and Herbivore density when controlling for year-year variation with a low R2
  • MODERATE significance in how the strength of the relationship between Spodoptera growth rate and Herbivore density varies among years with moderate R2

3.2- Year to year variation in Herbivore density and Spodoptera growth rate

There is an interactive effect of Spodoptera growth rate * Year on Herbivore density as shown in the ANOVA analysis above. To visualize this, we will construct a scatter plot of Herbivore density ~ Spodoptera growth rate data and layer on top of this least square means + standard errors for both Spodoptera growth rate and Herbivore density.

CONCLUSIONS

Spodoptera growth rate clearly decreases each year, however herbivore density does not seem to display year-year patterns. Very clear negative year-year effects on Spodoptera are present which suggestive of the plants are becoming more and more resistant to them. However, the same story doesn’t hold as true for herbivore density as it seems to be lower in 2023 than 2021, but the pattern is not nearly as apparent.

3.3 Individual year relationships for Herbivore density ~ Spodoptera growth rate

Residual plots

ANOVA TABLES

ANOVA results log Herbivore count ~ log Avg Spod mass across all species
Year P_value Marginal_R2
2021 0.49 0.006
2022 0.62 0.002
2023 0.02 0.052

CONCLUSIONS Only 2023 shows a moderately significant relationship between Spodoptera growth rate and Herbivore density when controlling for species, however with a low R2 of 0.05.

4- Between-species analysis

We previously examined Spodoptera ~ Herbivore relationships at a community level, treating species as a random effect, however this relationship may exhibit scale-dependency. As such, we will now look at Spodoptera ~ Herbivore density relationships between species across years. To do this, we will calculate species means for Spodotpera growth rate and Herbivore density across years and construct a model to evaluate.

4.1- Species means table

Species means for log Herbivore count, log Spodoptera Mass across years
Species Year Mean_log_abundance_herb Mean_log_spodmass_mg
1 Acmispon_glaber 2021 4.255376 3.8633378
2 Artemisia_californica 2021 4.400411 0.6024247
3 Diplacus_aurantiacus 2021 2.189648 0.1831020
4 Encelia_californica 2021 2.717340 0.7764338
5 Eriogonum_fasciculatum 2021 2.191305 3.6684611
6 Grindelia_camporum 2021 4.558662 1.6960144
7 Isocoma_menziesii 2021 4.489899 0.9674563
8 Malacothamnus_fasciculatus 2021 5.519561 1.3252105
9 Malacothrix_saxatilis 2021 2.342781 1.8846512
10 Mirabilis_laevis 2021 3.384199 1.2507455
11 Salvia_apiana 2021 4.757645 2.4386826
12 Salvia_mellifera 2021 4.356426 2.5597448
13 Sisyrinchium_bellum 2021 1.856765 0.4100945
15 Acmispon_glaber 2022 1.045045 2.1520206
16 Artemisia_californica 2022 3.758577 -0.4407271
17 Diplacus_aurantiacus 2022 3.242293 -0.9097240
18 Encelia_californica 2022 1.546300 -0.1885777
19 Eriogonum_fasciculatum 2022 3.856301 0.5842091
20 Grindelia_camporum 2022 2.225011 1.2479602
21 Isocoma_menziesii 2022 2.976124 -0.7357583
22 Malacothamnus_fasciculatus 2022 4.251484 0.6008157
23 Malacothrix_saxatilis 2022 1.754840 -0.3156419
24 Mirabilis_laevis 2022 2.218370 -0.6121919
25 Salvia_apiana 2022 4.467119 -0.3342514
26 Salvia_mellifera 2022 2.882546 0.0796477
27 Sisyrinchium_bellum 2022 1.118831 -1.2634918
29 Acmispon_glaber 2023 6.313361 NA
30 Artemisia_californica 2023 4.248330 0.1732868
31 Diplacus_aurantiacus 2023 2.916212 -2.3025851
32 Encelia_californica 2023 4.324319 -0.7675284
33 Eriogonum_fasciculatum 2023 2.842681 -1.1137470
34 Grindelia_camporum 2023 3.433987 1.6094379
35 Isocoma_menziesii 2023 2.968976 -2.3025851
36 Malacothamnus_fasciculatus 2023 3.812132 -0.9878210
37 Malacothrix_saxatilis 2023 1.832662 -0.8864200
38 Mirabilis_laevis 2023 2.419181 -2.3025851
39 Salvia_apiana 2023 2.935710 -0.6737646
40 Salvia_mellifera 2023 3.782973 -2.0110518
41 Sisyrinchium_bellum 2023 1.288061 -2.3025851

4.2- ANOVA of Herbivore density ~ Spodoptera growth rate

The model we will run is as follows. NOTE: The fixed effects are species means rather than individual values. Also, Species is included as a random effect due to repeated measures each year.

Herbivore density ~ Spodoptera growth rate + Year + (1|Species)

ANOVA results of Herbivore density ~ Spodoptera growth rate
Chisq Df Pr(>Chisq) Marginal_R2
Mean_log_spodmass_mg 0.1413087 1 0.7069834 0.0699179
Year 4.6282123 2 0.0988545 0.0949062

Plot of species means for Herbivore density ~ Spodoptera growth rate

TAKEAWAY

As seen from both the ANOVA and the scatter plot, there are no consistent patterns of Spodoptera ~ Herbivore density between species across years. For example, some show a positive relationship with similar slopes such as Artemisia californica, Malacothamnus fasciculatus and Salvia apiana whereas others show either a very steep positive slope (Grindelia camporum) or a negative slope (Diplacus aurantiacus, Eriogonum fasciculatum).

4.3- Year to year variation in Herbivore density ~ Spodoptera growth rate

Since there are strong year-year effects between species, we used the following model for analyses using subsets of the dataset for each year.

Herbivore density ~ Spodoptera growth rate + (1|Species)

2021

2021 ANOVA results of Herbivore density ~ Spodoptera growth rate
Chisq Df Pr(>Chisq) Marginal_R2
Mean_log_spodmass_mg 0.3970256 1 0.5286294 0.0320258739838356

2022

2022 ANOVA results of Herbivore density ~ Spodoptera growth rate
Chisq Df Pr(>Chisq) Marginal_R2
Mean_log_spodmass_mg 0.0350246 1 0.8515441 0.00268696329294603

2023

2023 ANOVA results of Herbivore density ~ Spodoptera growth rate
Chisq Df Pr(>Chisq) Marginal_R2
Mean_log_spodmass_mg 2.718746 1 0.0991759 0.198177469452273

5- Within-species analysis

We know the relationship between Spodoptera growth rate and Herbivore density does not show a signficiant relationship between species, but we dont yet know the strength of the relationship between Spodoptera growth rate and Herbivore density within a species. This next set of analyses will explore this in detail.

Analyses will be conducted at the species level and here focus on within-species variation across years with Plot and Plant as random effects.

The generalized linear mixed model (lmer) that we will use to fit the data is:

Herbivore density ~ Spodoptera growth rate + Species + Spodoptera growth rate * Species + Year + (1|Plot) + (1|Plant)

5.1- ANOVA of Herbivore count ~ Spodoptera growth rate

Similar to 3.1, we will look at the relationship between Spodoptera growth rate and Herbivore density but this time with Species as a fixed effect.

ANOVA results of log Herbivore count ~ log Avg Spod mass
Chisq Df Pr(>Chisq) Marginal_R2
log(Avg_spod_mass_mg) 0.4420116 1 0.5061531 0.0478649
Species 117.3201822 12 0.0000000 0.3191793
Year 29.7504088 2 0.0000003 0.4203546
log(Avg_spod_mass_mg):Species 39.4562448 12 0.0000885 0.0557591

RESULTS

  • NO significance in relationship between Spodoptera growth rate and Herbivore density when when controlling for year-year variation and species effects.
  • STRONG significance between species in Herbivore density
  • STRONG significance between years in Herbivore density, similar to what was seen in 3.1
  • STRONG significance in the interactive effect of Species and Spodoptera growth rate on Herbivore density.

5.2- Herbivore density ~ Spodoptera growth rate within Species

We will now explore the interactive effect by addressing how the effect of Spodoptera growth rate on Herbivore density varies across species. To do this, we will conduct a species by species analysis using data for all years then create a table to display the data. Model used is as follows:

Herbivore density ~ Spodoptera growth rate + Year + (1|Plot) + (1|Plant)

ANOVA Results

ANOVA table for Species-level effects of Spodoptera growth rate on Herbivore density across years
Species P_value Estimate_spodoptera Std_error R2
Acmispon_glaber 0.002 -0.620 0.200 0.07
Artemesia_californica 0.660 0.070 0.180 0.07
Diplacus_aurantiacus 0.880 -0.060 0.360 0.10
Encelia_californica 0.490 -0.100 0.150 0.05
Eriogonum_fasciculatum 0.004 0.330 0.120 0.03
Grindelia_camporum 0.240 0.670 0.570 0.25
Isocoma_menziesii 0.430 -0.140 0.170 0.71
Malacothamnus_fasciculatus 0.010 0.210 0.080 0.34
Malacothrix_saxatilis 0.050 0.840 0.430 0.18
Mirabilis_laevis 0.860 -0.040 0.210 0.24
Salvia_apiana 0.910 0.009 0.085 0.22
Salvia_mellifera 0.004 0.390 0.130 0.09
Sisyrinchium_bellum 0.370 0.360 0.410 0.03

5.3- Year to year variation in Herbivore density

Similar to 3.3, we will run a model and extract least squared means to look at year to year variation in herbivore density but this time with species as a fixed effect, rather than a random effect. This will allow us to understand the main effect of year on Herbivore density.

ANOVA results of Herbivore density ~ Year
Chisq Df Pr(>Chisq)
Year 45.72114 2 0
Least Square means of Herbivore density ~ Year
Year lsmean SE df lower.CL upper.CL
2021 3.579405 0.0945521 575 3.393695 3.765115
2022 2.761432 0.1030126 575 2.559105 2.963758
2023 3.143610 0.1103002 575 2.926969 3.360250

Nell, Colleen S., and Kailen A. Mooney. 2019. “Plant Structural Complexity Mediates Trade-Off in Direct and Indirect Plant Defense by Birds.” Ecology 100 (10). https://doi.org/10.1002/ecy.2853.
Pratt, Jessica D., Andrew Datu, Thi Tran, Daniel C. Sheng, and Kailen A. Mooney. 2017. “Genetically Based Latitudinal Clines in Artemisia Californica Drive Parallel Clines in Arthropod Communities.” Ecology 98 (1): 79–91. https://doi.org/10.1002/ecy.1620.
You, Yanrong, Chunpeng An, and Chuanyou Li. 2020. “Insect Feeding Assays with Spodoptera Exigua on Arabidopsis Thaliana.” BIO-PROTOCOL 10 (5). https://doi.org/10.21769/BioProtoc.3538.