DOI: 10.5061/dryad.5dv41nsh2
title: Nutrient addition, but not predator exclusion, shapes arthropod communities and herbivory in a temperate forest

description: Plants support diverse arthropod communities, and arthropod herbivores
 respond differently to plant traits, nutritional content, and defences,
 which influence their host plants selection, survival and performance
 (bottom-up control). At the same time, arthropod herbivores are affected
 by their interactions with predators (top-down control). Investigating how
 these forces interact, and how they are affected by nutrient availability,
 is crucial for understanding the mechanisms driving herbivore populations
 and their impact on ecosystems. We investigated how predator exclusion and
 fertilisation affect arthropod densities, sizes, and herbivory on two
 temperate forest tree species. Using a factorial design, we compared
 fertilised and unfertilised trees with and without the exclusion of flying
 vertebrate predators in the forest understory during September 2020 and
 2021. We collected and identified arthropods into feeding guilds.
 Fertilisation, but not predator exclusion, increased herbivory damage on
 trees as well as the size of predatory arthropods on the fertilised trees.
 Nutrient addition and predator exclusion had no significant effect on
 arthropod density. These patterns may indicate that the additional
 nutrients could have attracted herbivores, which in turn attracted their
 predators and may have enhanced their activity, thus potentially
 offsetting detectable changes in herbivore density. These results suggest
 that nutrient enrichment and predator exclusion interact, with nutrient
 addition affecting plant growth and herbivory damage primarily, but also
 increasing the size of predatory arthropods, but not the size of all
 arthropods. Effects of predator exclusion were less pronounced,
 potentially due to larger predatory arthropods compensating for the
 absence of flying vertebrate predators. Our study provide fully factorial
 field tests of top-down and bottom-up forces in a temperate forest
 understory, underscores the critical need to evaluate how diverse
 ecological interactions mediate the synergistic effects of nutrient
 pollution and the ongoing decline of insectivorous vertebrate predators on
 arthropod communities and herbivory damage, particularly as these
 preassures intensify in our rapidly changing environment.
description: <em><strong>Arthropod
 survey</strong></em> We used beating sheets
 to collect all arthropods from each of the trees. We measured all
 arthropods (to the nearest mm) and identified them into morphospecies,
 assigned them to an order or family if necessary for their assignment to
 one of the following four feeding guilds: predator, leaf chewer, leaf
 sucker, and as having “no relationship” (i.e., arthropod with no
 consumptive relationship with other arthropods or plants). Considering
 that the adult arthropod may belong to a different feeding guild than the
 juvenile, the developmental stage was considered in the guild
 assignment.  <em><strong>Herbivory damage
 survey</strong></em> Individual leaves with
 a piece of cm-scale were scanned (EPSON, 600 dpi, colour tiff format)
 within 12 hours after collection. To analyse herbivory, we first outlined
 the missing part of the leaves based on the expected shape in Photoshop®
 (following Sam et al., 2020). Then we calculated the remaining area (a)
 and the full expected area of each leaf (b, both in
 cm<sup>2</sup>) using ImageJ version 1.47 (National Institute
 of Health, USA), and we calculated the area eaten by the herbivores (c = b
 - a). We then calculated the proportion of leaf area lost as c/b for each
 leaf, which we used in our analyses as a mean across the 10 leaves per
 individual branch. For visualisation, we used percentages of herbivory.
 Next, we estimated the total number of leaves for each tree by taking the
 average of three independent estimates. We then determined the total leaf
 area per tree by multiplying the estimated number of leaves by the mean
 leaf area calculated from the leaves we used for the herbivory
 measurement. <em><strong>Leaf toughness
 survey</strong></em> We collected 10 random
 fresh leaves from each sapling. We measured leaf toughness using a FL50
 penetrometer with a Stepper Motor Powered Test Stand TVO 500N55S (Sauter
 GmbH) and a 1.5-mm-wide tip. We measured toughness once for each leaf,
 avoiding leaf veins (Onoda et al., 2011)
description: # **Nutrient addition, but not predator exclusion, shapes arthropod
 communities and herbivory in a temperate forest** Dataset DOI:
 [10.5061/dryad.5dv41nsh2](10.5061/dryad.5dv41nsh2) **Description of the
 data and file structure** This dataset was collected in a riparian forest
 (*Alnion incanae*) near Čejkovice, Czech Republic, as part of a study
 investigating the effects of nutrient addition and predator exclusion on
 arthropod communities and herbivory in a temperate forest. Using a fully
 factorial field experiment, we manipulated fertilization levels
 (fertilized, non-fertilized) and access to flying vertebrate predators
 (birds and bats) on two tree species: *Tilia cordata* and *Alnus
 glutinosa*. Over two years (2020–2021), we surveyed arthropods on
 experimental trees by collecting individuals and classifying them into
 morphospecies and feeding guilds (chewers, sap-suckers, predators, and
 arthropods with no direct relationship to plants or other arthropods). We
 measured arthropod density and body size, assessed herbivory damage by
 quantifying leaf area loss, and measured leaf toughness as a potential
 defense response. Our treatments had minimal to no impact on arthropod
 density. While we initially hypothesized that nutrient-rich foliage would
 attract more arthropods, our results did not support this assumption.
 Instead, we found that fertilized trees hosted larger herbivorous
 arthropods, which were subsequently followed by larger predatory
 arthropods. Additionally, fertilization—but not predator exclusion—led to
 increased herbivory damage, while trees where predators were removed
 exhibited lower leaf toughness compared to those with free predator
 access. **Files and variables** **File: arthropod_abundances.R**
 **Description:**\ This script conducts statistical analyses on arthropod
 abundance and density across different plant species and experimental
 treatments. It includes data preprocessing, model fitting using
 generalized linear mixed-effects models (GLMMs), model selection using
 Akaike Information Criterion (AICc), pairwise comparisons using estimated
 marginal means (EMMeans), and visualization of results. In addition, the
 script also performs Herbivore : Predator ratio analysis to quantify the
 balance between herbivore and predator pressure, and Guild composition
 analysis (PERMANOVA and stacked barplots) to test for community-level
 differences among treatments and species. **Variables:** * **code**:
 Identifier for each plant individual * **species**: Tree species where
 arthropods were collected (Tilia cordata, Alnus glutinosa) * **chewing**:
 Abundance of leaf-chewing arthropods * **sucking**: Abundance of
 sap-sucking arthropods * **pred**: Abundance of predatory arthropods *
 **no_relationship**: Arthropods with no clear relationship to plants or
 other arthropods * **not_sure**: Arthropods with uncertain feeding guild
 classification * **total_abundance**: Total arthropod abundance per plant
 individual * **treatment**: Experimental treatment (Control, No Predators,
 Fertilized, No Predators + Fertilized) * **average_LE**: Average leaf
 estimation per sample (calculated from two or three leaf estimations) *
 **average_leaf_size_cm2**: Mean size of individual leaves *
 **total_leaf_area_cm2**: Total leaf area of individual tree in cm^2^ *
 **total_leaf_area_m2**: Total leaf area of individual tree (in square
 meters) * **abundance_individuals_per_m2**: Standardized arthropod
 abundance per individual tree area (per m², density) * **year**: Year of
 data collection (2020, 2021) **File: abundances.xlsx** **Description:**\
 This dataset is used with "Arthropod_abundances.R" to conduct
 statistical analyses on arthropod abundance and density across different
 plant species and experimental treatments. **Table columns:** * **code**:
 Identifier for each plant individual * **species**: Tree species where
 arthropods were collected (Tilia cordata, Alnus glutinosa) * **chewing**:
 Abundance of leaf-chewing arthropods * **sucking**: Abundance of
 sap-sucking arthropods * **pred**: Abundance of predatory arthropods *
 **no_relationship**: Arthropods with no clear relationship to plants or
 other arthropods * **not_sure**: Arthropods with uncertain feeding guild
 classification * **total_abundance**: Total arthropod abundance per plant
 individual * **treatment**: Experimental treatment (Control, No Predators,
 Fertilized, No Predators + Fertilized) * **average_LE**: Average leaf
 estimation per sample (calculated from two or three leaf estimations) *
 **average_leaf_size_cm2**: Mean size of individual leaves in square
 centimeters * **total_leaf_area_cm2**: Total leaf area of individual tree
 in square centimeters * **total_leaf_area_m2**: Total leaf area of
 individual tree (in square meters) * **abundance_individuals_per_m2**:
 Standardized arthropod abundance per individual tree area (per m²,
 density) * **year**: Year of data collection (2020, 2021)   **File:
 body_sizes.R** **Description:**\ This script conducts statistical analyses
 on arthropod body sizes across different treatments, plant species, and
 feeding guilds. It includes data preprocessing, model fitting using
 generalized linear mixed-effects models (GLMMs), model selection using
 Akaike Information Criterion (AICc), pairwise comparisons using estimated
 marginal means (EMMeans), and visualization of results. **Variables:** *
 **treatment**: Control, No Predators, Fertilized, No Predators +
 Fertilized * **species**: Plant species (Tilia cordata, Alnus glutinosa) *
 **guild**: Feeding guilds of arthropods (Chewers, Predators, Sapsuckers,
 No Relationship) * **individual**: Unique identifier for each plant *
 **year**: Year of data collection (2020-2021) * **bodylength_mm**:
 Measured body length of arthropods (in milimeters) **Dataset:
 sizes.xlsx**\ Description: This dataset is used with
 "body_sizes.R" to conduct analyses on arthropod body sizes
 across different treatments, plant species, and feeding guilds. Contains
 one table. **Table columns:** * **treatment**: Experimental treatment
 (Control, No Predators, Fertilized, No Predators + Fertilized) *
 **species**: Plant species (Tilia cordata, Alnus glutinosa) *
 **individual**: Unique identifier for each plant * **notes id**: More
 detailed description of the ID if necessary * **order**: Arthropod order *
 **notes_id:** Another details important for the ID * **guild**: Feeding
 guild of the arthropod * **year**: Year of data collection (2020-2021) *
 **bodylength_mm**: Measured body length of arthropods (in milimeters) *
 **stage**: Life stage of the arthropod **File: herbivory.R**
 **Description:**\ This script conducts statistical analyses on herbivory
 damage across different plant species and experimental treatments. It
 includes data preprocessing, model fitting using generalized linear
 mixed-effects models (GLMMs) with a beta error distribution, model
 selection using Akaike Information Criterion (AICc), pairwise comparisons
 using estimated marginal means (EMMeans), and visualization of results. It
 uses the dataset **herbivory.xlsx**. **Variables:** * **Treatment**:
 Experimental treatment (Control, No Predators, Fertilized, No Predators +
 Fertilized) * **Species**: Tree species where herbivory was measured
 (Tilia cordata, Alnus glutinosa) * **Season**: Year of data collection
 (2020, 2021) * **IndividualCode**: Unique identifier for sampled
 individual plants within the plant species and treatments * **Twig**:
 Identifier for individual twigs sampled per individual plants (Red, White,
 Blue mark at each plant individual) * **Individual**: Identifier for each
 plant per individual treatments and plant species * **Area.x_cm2**: The
 whole leaf without missing parts (in square centimeters) * **Area.y_cm2**:
 The area of the leaf with missing parts (in square centimeters) *
 **Herbivory_cm2**: Area.x - Area.y. The resulting herbivory on individual
 leaves (in square centimeters) * **HerbProp**: Proportion of leaf area
 lost due to herbivory * **herbPerc_%**: Percentage of leaf area lost due
 to herbivory in percents **Dataset: herbivory.xlsx**\ Description: This
 dataset is used with "herbivory.R" for statistical analyses on
 herbivory damage across different plant species and experimental
 treatments. Contains one table. **Table columns:** * **Treatment**:
 Experimental treatment (Control, No Predators, Fertilized, No Predators +
 Fertilized) * **Species**: Tree species where herbivory was measured
 (Tilia cordata, Alnus glutinosa) * **Season**: Year of data collection
 (2020, 2021) * **IndividualCode**: Unique identifier for sampled
 individual plants within the plant species and treatments * **Twig**:
 Identifier for individual twigs sampled per individual plants (Red, White,
 Blue mark at each plant individual) * **Individual**: Identifier for each
 plant per individual treatments and plant species * **Area.x_cm2**: The
 area of the leaf without missing parts (in square centimeters) *
 **Area.y_cm2**: The area of the leaf with missing parts (in square
 centimeters) * **Herbivory_cm2**: Area.x - Area.y. The resulting herbivory
 on individual leaves in square centimeters * **HerbProp**: Proportion of
 leaf area lost due to herbivory * **herbPerc_%**: Percentage of leaf area
 lost due to herbivory in percents   **File: Leaf_toughness(1).R**
 **Description:**\ This script conducts statistical analyses on leaf
 toughness across different plant species and experimental treatments. It
 includes data preprocessing, model fitting using linear mixed-effects
 models (LMMs), model selection using Akaike Information Criterion (AICc),
 pairwise comparisons using estimated marginal means (EMMeans), and
 visualization of results. **Variables:** * **Treatment**: Experimental
 treatment (Control, No Predators, Fertilized, No Predators + Fertilized) *
 **Species**: Tree species where herbivory was measured (Tilia cordata,
 Alnus glutinosa) * **Year**: Year of data collection (2020, 2021) *
 **Individual**: Unique identifier for sampled individual plants within the
 plant species and treatments * **m1_N-m10_N**: Measured pressure needed to
 penetrate the leaf in newtons (N) **Dataset: Leaf_toughness.xlsx**\
 Description: This dataset is used with "Leaf_toughness.R" to
 analyze leaf toughness across different plant species and experimental
 treatments. It includes one table. **Table columns:** * **Treatment**:
 Experimental treatment (Control, No Predators, Fertilized, No Predators +
 Fertilized) * **Species**: Tree species where herbivory was measured
 (Tilia cordata, Alnus glutinosa) * **Year**: Year of data collection
 (2020, 2021) * **Individual**: Unique identifier for sampled individual
 plants within the plant species and treatments * **m1_N-m10_N**: The
 measured pressure needed to penetrate the leaf in newtons (N)  
 **Code/software** **Arthropod_abundances.R** This script conducts
 statistical analyses on **arthropod density** across different plant
 species and experimental treatments. It includes data preprocessing, model
 fitting using **generalized linear mixed-effects models (GLMMs)**, model
 selection using **Akaike Information Criterion (AICc)**, pairwise
 comparisons using **estimated marginal means (EMMeans)**, and
 visualization of results. It also performs:\ • Herbivore : Predator ratio
 analysis to quantify the balance between herbivore and predator pressure,
 and\ • Guild composition analysis (PERMANOVA and stacked barplots) to test
 for community-level differences among treatments and species.\ It uses the
 dataset **“abundances.xlsx.”**It uses the dataset
 "abundances.xlsx" **Language and Environment** * **R Environment
 for Statistical Computing** * **Version:** R 4.3.1 **Dependencies** This
 analysis requires the following R packages: * **readxl** – For importing
 Excel data * **dplyr** – For data manipulation * **ggplot2** – For data
 visualization * **ggpubr** – For arranging multiple plots * **glmmTMB** –
 For running generalized linear mixed-effects models (GLMMs) * **lme4** –
 For fitting linear mixed-effects models * **bbmle** – For model selection
 using AICc * **emmeans** – For pairwise comparisons of estimated marginal
 means * **logspline** – For distribution fitting * **fitdistrplus** – For
 distribution fitting * **lmtest** – For likelihood ratio tests * **vegan**
 – For multivariate analyses (PERMANOVA using Bray–Curtis distances) *
 **reshape2** – For reshaping data into long format for stacked barplots  
 **Steps Performed in the Script** The following steps (1–8) describe
 analyses ranging from overall arthropod abundance to guild-level patterns
 and community composition. **1. Import and Prepare Data** * Loads the
 **abundances.xlsx** dataset containing arthropod abundance measurements
 across different treatments and plant species. * Converts necessary
 variables to factors: * **Treatment** (Control, No Predators, Fertilized,
 No Predators + Fertilized) * **Species** (*Tilia cordata*, *Alnus
 glutinosa*) * **Year** (Year of data collection) * **Code** (Unique
 identifier for each sample) * Converts **abundance** to numeric format and
 adjusts formatting inconsistencies. * Ensures factor variables are
 properly structured for analysis. **2. Fit Generalized Linear
 Mixed-Effects Models (GLMMs)** * Uses the **glmmTMB** package to model
 **arthropod density** as a function of treatment and species while
 controlling for year and sample code as random effects. * Runs multiple
 models: * **Null model** (M1.null) – Baseline model with no predictors. *
 **Treatment model** (M1.Treatment) – Tests the effect of fertilization and
 predator exclusion. * **Species model** (M1.species) – Tests the effect of
 plant species. * **Interaction model** (M1.Treatment.species_interaction)
 – Tests interaction between species and treatment. * **Full model**
 (M1.full) – Includes all predictors. **3. Model Selection** * Uses
 **AICctab()** from the **bbmle** package to compare model fits and select
 the best-performing model based on AICc values. * Performs **likelihood
 ratio tests (lrtest)** to determine statistical significance between
 models. **4. Estimated Marginal Means (EMMeans) Analysis** * Calculates
 **pairwise comparisons** of treatment and species effects using the
 **emmeans** package. * Extracts contrasts and response estimates for
 further analysis. * Saves results as CSV files. **5. Visualization** *
 Creates **scatter plots and error bar plots** using **ggplot2** to
 display: * Individual data points for arthropod abundance. * Estimated
 marginal means with confidence intervals. * Treatment-based color coding.
 * Saves the final plot as a high-resolution image. **6. Analysis for
 Feeding Guilds (Chewers and Predators)** **Analysis for Feeding Guilds**
 **Chewers (Leaf-Chewing Arthropods)** * Computes chewer density per m² of
 foliage, in the previous analysis is this already done in the data: *
 arth$abundanceCHEW &lt;- (arth$chewing / arth$total_leaf_area_m) + 0.1 *
 Fits GLMMs using a Gaussian family with log link: * **M2.null:** Baseline
 model with year as a fixed effect and sample code as a random effect. *
 **M2.Treatment:** Tests treatment effects. * **M2.species:** Tests species
 effects. * **M2.Treatment.species_interaction:** Tests interaction
 effects. * **M2.full:** Includes all predictors. * Performs model
 selection using AICctab(). * Uses emmeans to compute pairwise species
 comparisons and saves results. * Creates visualizations using ggplot2,
 showing: * Data points for chewer density. * Estimated marginal means with
 confidence intervals. **Predators (Predatory Arthropods)** * Computes
 predator density per m² of foliage, simirarly as for chewers: *
 arth$abundancePRED &lt;- (arth$pred / arth$total_leaf_area_m) + 0.000001 *
 Fits GLMMs using a Gaussian family with log link: * **M3.null:** Baseline
 model. * **M3.Treatment:** Tests treatment effects. * **M3.species:**
 Tests species effects. * **M3.Treatment.species_interaction:** Tests
 interaction effects. * **M3.full:** Includes all predictors. * Performs
 model selection using AICctab() and likelihood ratio tests. * Uses emmeans
 for pairwise species comparisons and saves results * Creates
 visualizations using ggplot2, showing: * Data points for predator density.
 * Estimated marginal means with confidence intervals. **Combined Plot for
 Chewers and Predators** * Arranges chewer and predator density plots side
 by side using ggpubr. * Saves the final combined plot. **7. Herbivore :
 Predator Ratio** •  Computes the herbivore : predator ratio per individual
 tree to assess the balance between herbivore and predator pressure. •  The
 ratio is calculated as (chewers / predators) standardized by total leaf
 area and then log-transformed for normality. •  Fits a GLMM (Gaussian
 family) testing for species-level differences in this ratio across years.
 •  Estimated marginal means (EMMeans) are computed for species effects. • 
 Results are visualized using ggplot2, showing both raw data points and
 estimated means with confidence intervals. •  The final figure is saved as
 “herbivore_predator_ratio_species.jpeg”.   **8. Guild Composition Analysis
 (PERMANOVA and Stacked Barplot)** •  Calculates relative abundance (%) of
 each arthropod feeding guild (chewers, suckers, predators,
 no_relationship) across treatments and tree species. •  Uses multivariate
 PERMANOVA (adonis2, Bray–Curtis distances) to test for differences in
 guild composition among treatments and species. •  Data are visualized as
 stacked barplots showing the proportion of each guild by treatment for
 both species. •  Adds text labels to display mean relative abundances per
 guild within treatments. •  The final figure is saved as
 “guild_composition_treatments.jpeg”. **Leaf_toughness.R** this script
 conducts statistical analyses on leaf toughness across different plant
 species and experimental treatments. It includes data preprocessing, model
 fitting using linear mixed-effects models (LMMs), model selection using
 Akaike Information Criterion (AICc), pairwise comparisons using estimated
 marginal means (EMMeans), and visualization of results. It uses the
 dataset *Leaf_toughness.xlsx*. **Language and Environment** **R
 Environment for Statistical Computing**\ Version: R 4.3.1 **Dependencies**
 This analysis requires the following R packages: * **readxl** – For
 importing Excel data * **dplyr** – For data manipulation * **reshape2** –
 For reshaping to long data version * **ggplot2** – For data visualization
 * **bbmle** – For model selection using AICc * **lme4** – For fitting
 linear mixed-effects models * **emmeans** – For pairwise comparisons of
 estimated marginal means   **Steps Performed in the Script** **1. Import
 and Prepare Data** * Loads the *Leaf_toughness.xlsx* dataset containing
 leaf toughness measurements across different treatments and plant species.
 * Reshapes the dataset into a long format using melt() from the reshape2
 package. * Converts necessary variables to factors: * **Treatment**
 (Control, Fertilized, etc.) * **Species** (*Tilia cordata, Alnus
 glutinosa*) * **Year** (Year of data collection) * **Individual** (Unique
 identifier for each sample) * Converts response variables to numeric
 format and ensures proper structuring for analysis. * Removes missing
 values and zeros from the dataset. * filters out *Quercus* species which
 we did not use later in the analyses. **2. Fit Linear Mixed-Effects Models
 (LMMs)** * Uses the lmer() function from the **lme4** package to model
 leaf toughness as a function of treatment and species while controlling
 for individual variation as a random effect. * Runs multiple models: *
 **Null model (mixed_model.0)** – Baseline model with no predictors. *
 **Treatment model (mixed_model.T)** – Tests the effect of fertilization
 and other treatments. * **Species model (mixed_model.S)** – Tests the
 effect of plant species. * **Additive model (mixed_model.TS)** – Includes
 species and treatment as fixed effects. * **Interaction model
 (mixed_model.TiS)** – Tests the interaction between species and treatment.
 **3. Model Selection** * Uses AICtab() from the **bbmle** package to
 compare model fits and select the best-performing model based on AICc
 values. * Performs likelihood ratio tests (lrtest()) to determine
 statistical significance between models. **4. Estimated Marginal Means
 (EMMeans) Analysis** * Uses the emmeans() function to compute pairwise
 comparisons of treatment effects within each species. * Saves results as
 CSV files. **5. Visualization** * Creates scatter plots and error bar
 plots using **ggplot2** to display: * Individual data points for leaf
 toughness measurements. * Estimated marginal means with confidence
 intervals. * Treatment-based color coding. * Saves the final plot as a
 high-resolution image.   **Herbivory.R** This script analyzes herbivory
 damage across different species, treatments, and seasons using generalized
 linear mixed models (GLMMs). It includes data preprocessing, model
 selection using AICc, estimated marginal means (EMMeans) calculations, and
 visualization of treatment effects. The dataset *herbivory.xlsx* is used
 for analysis. **Language and Environment** **R Environment for Statistical
 Computing**\ Version: R 4.3.1 **Dependencies** The script requires the
 following R packages: * **readxl** – For importing Excel data * **dplyr**
 – For data manipulation * **ggplot2** – For data visualization *
 **glmmTMB** – For fitting generalized linear mixed models * **bbmle** –
 For model selection using AICc * **emmeans** – For pairwise comparisons of
 estimated marginal means * **lmtest** – For likelihood ratio tests *
 **ggpubr** – For combining plots   **Steps Performed in the Script** **1.
 Import and Prepare Data** * Loads the *herbivory.xlsx* dataset. * Converts
 key categorical variables to factors: * **Treatment** (control, no
 predators, fertilized, no predators + fertilized) * **Species** (Plant
 species) * **Season** (Season of the experiment, 2020, 2021) *
 **Individual** (Unique identifier for each plant) * **Twig** (Leaf cluster
 identifier) * Filters the dataset to include only the four main treatment
 groups. * Computes summary statistics for herbivory percentage and
 calculates means, standard deviations, and standard errors. **2. Fit
 Generalized Linear Mixed Models (GLMMs)** * Uses glmmTMB() to fit beta
 regression models with a logit link function. * Models herbivory
 proportion as a function of treatment, species, and season, including
 random effects for individual plants. * Fits multiple models: * **Null
 model (H1.null)** – Baseline model with no predictors. * **Treatment model
 (H1.Treatment)** – Tests the effect of fertilization and predator
 exclusion. * **Species model (H1.Species)** – Tests species differences in
 herbivory. * **Interaction model (H1.Treatment.Species_interaction)** –
 Includes an interaction term between treatment and species. * **Full model
 (H3.full)** – Incorporates treatment, species, and season. **3. Model
 Selection** * Compares models using AICc with AICctab(). * Performs
 likelihood ratio tests (lrtest()) to evaluate model significance. *
 Selects the best-fitting model (H3.full). **4. Estimated Marginal Means
 (EMMeans) Analysis** * Uses emmeans() to estimate pairwise differences in
 herbivory among treatments within each species. **5. Visualization** *
 Creates a scatter plot with error bars using **ggplot2**: * Herbivory
 percentage is plotted by treatment and species. * Individual data points
 are jittered for readability. * Estimated marginal means are displayed
 with 95% confidence intervals. * Saves the final plot as a high-resolution
 image. * **Body_Sizes.R** This script analyzes insect body size variation
 across different guilds, species, and treatments. It fits generalized
 linear mixed models (GLMMs), selects the best model using AICc, computes
 estimated marginal means (EMMeans), and visualizes results. The dataset
 *sizes.xlsx* is used for analysis. **Language and Environment** **R
 Environment for Statistical Computing**\ Version: R 4.3.1 **Dependencies**
 The script requires the following R packages: * **readxl** – For importing
 Excel data * **dplyr** – For data preprocessing * **glmmTMB** – For
 fitting generalized linear mixed models * **bbmle** – For model selection
 * **emmeans** – For estimated marginal means analysis * **ggplot2** – For
 data visualization   **Steps Performed in the Script** **1. Import and
 Prepare Data** * Loads the *sizes.xlsx* dataset. * Converts key
 categorical variables to factors: * **Treatment** (control, no predators,
 fertilized, no predators + fertilized) * **Species** (Plant species) *
 **Guild** (Functional feeding group) * **Year** (Year of the sample
 collection) * **Individual** (Unique identifier for each insect) *
 **Order** (Taxonomic order of individual arthropods) * Removes missing
 values and unnecessary columns. **2. Fit Generalized Linear Mixed Models
 (GLMMs)** * Uses glmmTMB() to fit Gamma regression models with a log link
 function. * Models insect body size as a function of treatment, species,
 and guild, with year and individual as random effects. * Fits multiple
 models: * **Null model (sizes.null)** – Baseline model with no predictors.
 * **Treatment model (sizes.treat)** – Examines treatment effects. *
 **Guild model (sizes.guild)** – Examines guild-based differences. *
 **Species model (sizes.spec)** – Examines species-based differences. *
 **Full model (sizes.treat.guild.spec)** – Includes treatment, guild, and
 species as fixed effects. * **Interaction models** – Test interactions
 between treatment, guild, and species. **3. Model Selection** * Uses
 AICctab() to compare model fits based on AICc values. * Performs
 likelihood ratio tests (lrtest()) to determine statistical significance of
 model improvements. * Selects the best-performing model
 (sizes.treat.guild_interact). **4. Estimated Marginal Means (EMMeans)
 Analysis** * Uses emmeans() to estimate body size differences across
 guilds and treatments. * Saves results as CSV files. **5. Visualization**
 * Creates a scatter plot with error bars using **ggplot2**: * Individual
 data points for body size are plotted by guild and treatment. * Estimated
 marginal means with 95% confidence intervals are displayed. * Colors
 differentiate treatment groups. * Saves the final plot as a
 high-resolution image. **Missing or Unavailable Data** Cells containing
 'NA' indicate missing or non-applicable values as follows: * n/a
 – Data not available for that observation (e.g., measurement not taken,
 specimen lost, or unreadable sample). * 0 values – Represent true zero
 counts (e.g., no arthropods recorded).  These missing or unavailable data
 are not replaced or infilled, to avoid interfering with R scripts that
 handle them as NA or missing. All analyses explicitly include na.rm = TRUE
 or use complete.cases() to manage missing data. **Access information** *
  n/a