DOI: 10.5061/dryad.brv15dvm3
title: Regional variation in insularity effects on acorn herbivory in European oaks challenges predictions of lower herbivory on islands

description: Ecological theory predicts that herbivore abundance and diversity should
 be lower on islands compared to the mainland, resulting in lower
 herbivory. However, island-mainland comparisons have focused on herbivory
 in vegetative tissues, primarily leaves, while herbivory in reproductive
 tissues, such as fruits and seeds, which directly influence seedling
 establishment and plant fitness, has lagged behind. In this study, we
 compared insect herbivory in oak (Quercus) reproductive structures (i.e.,
 acorn damage) across seven species in three regions: Lesbos Island vs.
 mainland Greece, the Balearic Islands vs. mainland Spain, and Bornholm
 Island vs. mainland Sweden. For each species, we selected three islands
 and three mainland populations, and sampled 10 acorns from four trees per
 population (N = 150 trees). We also analysed acorn chemical defences
 (phenolics) and nutrients (nitrogen and phosphorus) to test their
 association with herbivory and their contribution to island-mainland
 differences in acorn damage. We obtained climatic data for each population
 to assess whether these factors mediated insularity effects on acorn
 traits and herbivory. We found no main effect of insularity, meaning acorn
 damage did not differ overall between mainland and island populations.
 However, we observed region-specific patterns whereby acorn damage was
 higher in the Balearic Islands than in mainland Spain, while no
 significant mainland-island differences were found in the other regions.
 There were no main insularity effects on acorn traits, but a
 region-by-insularity interaction affected phosphorus, with higher average
 values in mainland Sweden than in Bornholm Island, whereas the reverse was
 observed in Greece. A follow-up mechanistic test indicated that insularity
 indirectly influenced acorn traits (e.g., phosphorus) via climate, but
 these trait differences did not explain acorn damage patterns. Our results
 challenge the theory on insularity effects on herbivory and suggest that
 local processes are driving contrasting outcomes for acorn herbivory
 across the studied regions.
description: <em><strong>Field sampling and acorn damage
 assessment</strong></em> In September 2023,
 at the end of the growing season, we selected populations of oak species
 from both island and mainland locations. For each species, we chose three
 islands and three mainland populations, except for <em>Q.
 ithaburensis</em> and <em>Q. pubescens</em> in mainland
 Greece, for which we sampled one and two populations, respectively, as
 well as for <em>Q. suber</em> in the Balearic Islands, for
 which we sampled two populations. This resulted in a total of 38 oak
 populations. Each population consisted of at least 15 reproductive trees
 and was at least 2 km apart from any other selected population of the same
 species. From each population, we randomly selected four individual trees
 except for <em>Q. suber</em> in the Balearic Islands, for
 which we sampled two trees. The total number of sampled trees was 150.
 Then, for each tree, we haphazardly collected 10 acorns found below the
 tree canopy. The acorns were oven-dried at 40ºC for 48 hours. After
 drying, we assessed the presence of insect herbivore (e.g., weevils,
 moths) damage in each acorn (Xiao et al. 2007). Acorns were categorized as
 either "damaged” or “undamaged” based on the presence of larvae or
 exit holes (Xiao et al. 2007). We then calculated the proportion of acorns
 damaged by insect herbivores in each tree (“acorn damage” hereafter) as
 the ratio of damaged acorns to the total number of acorns, which was used
 in subsequent statistical analyses (see “Statistical analyses” section
 below). <em><strong>Acorn
 traits</strong></em> We analysed acorn
 traits potentially associated with acorn damage, specifically phosphorus
 and nitrogen content as proxies of acorn nutritional value, as well as the
 concentration of total phenolics as a proxy of chemical defences. We
 analysed two to three acorns per tree with little to no insect herbivore
 damage to avoid biases in trait levels due to local induction from insect
 damage (albeit effects of systemic induction on undamaged acorns cannot be
 discarded). To quantify nitrogen and phosphorous
 content, we ground together acorns from each tree using liquid nitrogen
 and digested 0.1 g of the ground acorn material in a mixture of
 selenic-sulfuric acid and hydrogen peroxide (Moreira et al., 2012).
 Diluted aliquots of the digestion were then analysed by colorimetry. We
 quantified nitrogen content using the indophenol blue method whereas
 phosphorus content was quantified with the molybdenum blue method, in both
 cases using a Biorad 650 microplate reader (Bio-Rad Laboratories,
 Philadelphia, PA, USA) at 650 nm and 700 nm, respectively (Walinga et al.,
 1995). We expressed nitrogen and phosphorus concentrations in mg g^-1^
 tissue on a dry weight basis. To quantify the
 concentration of total phenolic compounds, we used the same ground acorn
 material than for nutritional trait analyses. Briefly, we incubated 20 mg
 with 1 ml of 70% methanol in an ultrasonic bath for 15 minutes (Moreira et
 al., 2012). Total phenolic content was determined colorimetrically using
 the Folin-Ciocalteu assay, with a Biorad 650 microplate reader (Bio-Rad
 Laboratories, Philadelphia, PA, USA) at 740 nm. We quantified total
 phenolic content with a tannic acid standard curve and expressed the
 results as mg of tannic acid equivalents per gram of dry plant tissue (mg
 g^-1^) (Sampedro et al., 2011).
 <em><strong>Climate
 data</strong></em> We characterized the
 climatic conditions at each population using a subset of eight bioclimatic
 variables from the WorldClim database version 2.1 database at a 2.5-
 minute resolution (Fick and Hijmans., 2017). These variables were: BIO1
 (annual mean temperature, °C), BIO4 (temperature seasonality, expressed as
 the standard deviation of temperature among months × 100), BIO5 (maximum
 temperature of the warmest month, °C), BIO6 (minimum temperature of the
 coldest month, °C), BIO12 (annual precipitation, mm), BIO13 (precipitation
 of the wettest month, mm), BIO14 (precipitation of the driest month, mm),
 BIO15 (precipitation seasonality, expressed as standard deviation of
 precipitation across months). We summarized the population-level climatic
 data using principal component analysis (PCA) with the
 <em>FactorMineR</em> package  in R software (Husson et al.,
 2006). We retained the first principal component (PC1 climate, explaining
 61.28% of the climatic variation), which was positively correlated with
 temperature and precipitation seasonality, and negatively correlated with
 precipitation during the driest month (Fig. S2). We used the z-values from
 PC1 climate in our statistical analyses (see below).
 <em><strong>Statistical
 analyses</strong></em> We conducted all
 analyses using R version 4.2.1. First, we ran a generalized linear mixed
 model with a beta-binomial distribution (logit-link function) to test the
 effects of insularity (two levels: mainland or island), region (three
 levels: Lesbos Island-mainland Greece, Balearic Islands-mainland Spain,
 and Bornholm Island-mainland Sweden), and their interaction (all treated
 as fixed factors) on the proportion of damaged acorns per tree using the
 <em>glmmTMB</em> function from the
 <em>glmmTMB</em> package. We also included oak species and
 population within species as random factors. We assessed the significance
 of fixed effects using Wald χ² tests implemented via the
 <em>Anova</em> function from the <em>car</em>
 package and evaluated the contribution of random effects to the variance
 explained by comparing models with and without each random term using
 likelihood ratio tests (LRT) implemented via the
 <em>anova</em> function from the base
 <em>stats</em> package. In the event of a significant
 interaction (<em>p</em> &lt; 0.05) between insularity and
 region, we ran post-hoc mean contrasts for island vs. mainland within each
 region, using the <em>lsmeans</em> function from the
 <em>lsmeans</em> package, with Sidak correction for multiple
 comparisons and degrees of freedom estimated using the Kenward-Roger
 approximation. Pseudo-R² values for this model were calculated using the
 <em>performance</em> package. Second, we
 ran generalized linear mixed models with Gaussian error distribution to
 test the effects of insularity, region, and their interaction on each
 acorn trait, again using the <em>glmmTMB</em> function from
 the <em>glmmTMB</em> package. We also included species and
 population within species as random factors. We log-transformed phosphorus
 and nitrogen content to achieve normality of residuals. For these Gaussian
 models, we specified a variance structure to account for
 heteroscedasticity between insularity groups. Including this variance
 structure allowed to explicitly model group-specific residual variance,
 thereby improving model fit and ensuring valid estimation of fixed effects
 and standard errors. We applied the same procedure described above to
 assess fixed and random effects, using Wald χ² tests for fixed effects,
 and LRT for random effects. Again, in the case of significant interactions
 (<em>p</em> &lt; 0.05) between insularity and region, we
 ran post-hoc contrasts between means for island vs. mainland within each
 region, using the same <em>lsmeans</em> procedure as described
 above. For all models, we conducted residual diagnostics using the
 <em>DHARMa</em> package, and present these diagnostics in
 Figure S3. Next, we ran a piecewise structural equation
 model (pSEM) using population means for all variables in which we tested
 for associations between insularity, climatic variation (PC1 climate
 scores), acorn traits, and acorn damage. Specifically, we included the
 following direct effects: (1) of insularity, coded as a dummy variable
 (mainland = 0, islands = 1) on PC1 climate, acorn traits and acorn damage;
 (2) of PC1 climate on each acorn trait (nitrogen, phosphorous, and total
 phenolics) and on acorn damage (i.e., the latter to test for direct
 abiotic forcing on herbivores, an alternative to plant trait-mediated
 abiotic effects); and (3) of each acorn trait on acorn damage.
 Additionally, we also tested the indirect effects of (4) insularity on
 acorn traits, mediated by climate, and (5) insularity on acorn damage,
 mediated by acorn traits. The pSEM was implemented using the
 <em>psem</em> function from the
 <em>piecewiseSEM</em> package. The component models in the
 pSEM included Gaussian models and one non-Gaussian (binomial) model for
 acorn damage fitted with the <em>lmer</em> and
 <em>glmer</em> functions respectively, both from the
 <em>lme4</em> package. Direct effects were calculated as
 standardized path coefficients between variables, while indirect effects
 were calculated as the product of the intervening direct effects for the
 specified causal pathway. To estimate indirect effects with 95% confidence
 intervals, we used the <em>semEff</em> function from the
 <em>semEff</em> package, which applies bootstrap resampling to
 provide robust estimates.
description: # Acorn insect damage and associated traits across Island-Mainland regions
 [https://doi.org/10.5061/dryad.brv15dvm3](https://doi.org/10.5061/dryad.brv15dvm3) ## Description of the data and file structure We provide data on acorn damage by insect herbivores, acorn phenolic content, and acorn nutrient content, for oak (*Quercus*) acorns collected across 38 populations belonging to seven oak species distributed across three island regions and their mainland counterparts. We also provide data on climate conditions at each population.  ### Files and variables #### File: Coordenadas.csv **Description:** Coordinates for each study population ##### Variables * POP_ID: Population unique code * LAT: latitude (decimal) * LONG: Longitude (decimal) #### File: data.csv **Description:** Acorn damage data at the tree level ##### Variables * ID: tree unique ID * INS: Insularity (whether island or continent) * SYST: Study region * SP: Study oak species * POP: population nested within species and region  * DAMAGE: number of damage acorns * NO DAMAGE: number of undamaged acorns * Total: proportion of damaged acorns from the total collected #### File: data_nut.csv **Description:** Data on the concentration of phenolics, nitrogen, and phosphorus in acorn content at the tree level ##### Variables * ID: tree unique ID  * INS: Insularity (whether island or continent) * SYST: Study region * SP: Study oak species * POP: population nested within species and region  * WEIGHT: weighed dry material for analyses (g) * mgP/gleaf: concentration of phosphorus (mg/ g) * mgN/gleaf: concentration of nitrogen (mg/ g) * Phen: concentration of phenolics (mg/ g) #### File: climate_data.csv **Description:** Climate data for the 8 bioclimatic variables extracted from World CLim ##### Variables * POP_ID: Unique population ID * wc2.1_2.5m_bio_1: annual mean temperature, °C * wc2.1_2.5m_bio_4: temperature seasonality, expressed as the standard deviation of temperature among months × 100 * wc2.1_2.5m_bio_5: maximum temperature of the warmest month, °C * wc2.1_2.5m_bio_6: minimum temperature of the coldest month, °C * wc2.1_2.5m_bio_12: annual precipitation, mm * wc2.1_2.5m_bio_13: precipitation of the wettest month, mm * wc2.1_2.5m_bio_14: precipitation of the driest month, mm * wc2.1_2.5m_bio_15:precipitation seasonality, expressed as standard deviation of precipitation across months ## Code/software The provided analysis was conducted using **R** and R Markdown. To reproduce the results, you need: * **R (version 4.x or later)** – The core statistical software used for data analysis. * **RStudio (optional, but recommended)** – A graphical interface for R that facilitates script execution and visualization. * **R Markdown** – The document format used for combining text, code, and output into a single report. #### Loaded Packages The analysis relies on the following R packages, which must be installed and loaded: * `lmerTest` – Linear mixed-effects models with p-values. * `lsmeans` – Least-squares means (for post hoc analysis). * `ggplot2` – Data visualization. * `tidyr` – Data manipulation. * `tidyverse` – Collection of data science packages. * `ggpubr` – Publication-ready plots. * `ggrepel` – Improved label placement for ggplot2. * `raster` – Geographic data manipulation. * `FactoMineR` – Multivariate analysis, including PCA. * `glmmTMB` – Generalized linear mixed models. * `DHARMa` – Residual diagnostics for mixed models. #### Workflow Description The workflow follows these steps: 1. **Loading and processing Data**: * Reads in **data.csv** (contains the main dataset). * Reads in **coordenadas.csv** (spatial coordinates of populations). * Merges coordinate data with the main dataset. 2. **Statistical Analysis & Visualization**: * Performs statistical modeling (likely mixed-effects models). * Generates plots (maps, data visualizations). ## Access information NA