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.

Field sampling and acorn damage assessment

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 Q. ithaburensis and Q. pubescens in mainland Greece, for which we sampled one and two populations, respectively, as well as for Q. suber 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 Q. suber 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).

Acorn traits

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).

Climate data

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 FactorMineR 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).

Statistical analyses

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 glmmTMB function from the glmmTMB 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 Anova function from the car 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 anova function from the base stats package. In the event of a significant interaction (p < 0.05) between insularity and region, we ran post-hoc mean contrasts for island vs. mainland within each region, using the lsmeans function from the lsmeans 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 performance 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 glmmTMB function from the glmmTMB 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 (p < 0.05) between insularity and region, we ran post-hoc contrasts between means for island vs. mainland within each region, using the same lsmeans procedure as described above. For all models, we conducted residual diagnostics using the DHARMa 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 psem function from the piecewiseSEM package. The component models in the pSEM included Gaussian models and one non-Gaussian (binomial) model for acorn damage fitted with the lmer and glmer functions respectively, both from the lme4 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 semEff function from the semEff package, which applies bootstrap resampling to provide robust estimates.