Nutrient addition, but not predator exclusion, shapes arthropod communities and herbivory in a temperate forest
Data files
Oct 20, 2025 version files 432.59 KB
-
abundances.xlsx
20.85 KB
-
arthropod_abundances.R
19.58 KB
-
body_sizes.R
7.87 KB
-
herbivory.R
7.56 KB
-
herbivory.xlsx
242.66 KB
-
Leaf_toughness.R
5.85 KB
-
Leaf_toughness.xlsx
15.83 KB
-
README.md
27.17 KB
-
sizes.xlsx
85.21 KB
Abstract
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.
Dataset DOI: 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 cm2
- 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 <- (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 <- (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
Arthropod survey
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.
Herbivory damage survey
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 cm2) 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.
Leaf toughness survey
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)