Data analyses
We used the interpolation procedure proposed by Husson & Josse (2016)
to impute the missing soil and climate data for three edaphic islands
out of 20. Soil properties were aggregated into four predictors: 1) mean
depth, 2) depth coefficient of variation (CV), 3) sandiness index, and
4) fertility. Mean soil depth and CV were calculated at the island level
using all the measurements taken at each island (i.e. five soil depth
measurements/sampled individual = ~125
measurements/island, total >2500 measurements). The
sandiness index was calculated from textural parameters as ∑sand /
∑(silt + clay), i.e. a higher proportion of sand produced higher values
of this index, implying less water retention capacity. Soil fertility
was identified as the first axis of a Principal Component Analysis (PCA;
first axis explaining ~49% of the variance, with
positive scores associated with more fertile soils; Supplemental
Material 1) ran on the eight soil nutrient parameters. Regarding
climate, we conducted a PCA on the 18 variables, and we selected the
first two axes to be used as predictors in the models. Specifically, the
first PCA axis was positively related to the annual mean of air and soil
temperature and negatively associated with the seasonality of those
variables (PC1clim; explaining ~29% of the variance).
The second PCA axis was positively related to annual mean soil moisture
and less extreme soil temperatures (PC2clim; explaining
~24% of the variance; Supplemental Material 1). For
insularity, we included the target effect. Trait values were averaged
(using the three values/species/island) for each species at the edaphic
island scale.
Before running the models, we controlled for collinearity among the
seven soil, climate and insularity variables through Spearman’s rho
correlation coefficient (with Bonferroni correction), and no
collinearity issues were detected (Supplemental Material 1). Model-wise,
we first ran linear mixed-effect models (LMMs) for clonal and non-clonal
species separately. In the LMMs, we set the trait average at the island
scale as the response variable, the seven environmental and insularity
variables as predictors (i.e. fixed effects; scaled and centered),
treating species identity as a random effect (informing on the magnitude
of species-specific responses). The variance explained by fixed effects
alone (marginal R2), and by fixed and random effect
together (conditional R2) was calculated. Because we
aimed at gaining insights into intraspecific responses to different
environmental and insularity conditions (and because this study had an
inherent explorative component), we examined plant trait-environment
links at single trait vs. single predictor level, using bivariate
ordinary least squares (OLS) linear models, setting species identity as
a grouping factor (i.e. interacting with the single predictor). We
identified important relationships based on the model coefficient and
its 95% confidence intervals (direction and robustness of the
relationship), R2 (goodness of fit and strength of the
relationship), p-value (significance of the relationship). Except for
storage tissue, all the other traits needed log-transformation to
accommodate the normality of data distribution and homoscedasticity of
model residuals, whereas environmental predictors did not require any
transformation. For clonal and non-clonal species separately, we
visualized the multivariate trait space identified by functional traits
and the species occupancy in this space (i.e. the non-acquisitive
persistence niche) through Nonmetric Multidimensional Scaling (NMDS). As
traits were measured in different units, we used Gower distance in the
NMDS and set 100 random starts and two dimensions. We conducted all the
analyses in R version 4.0.1 (R Core Team 2020) using functions available
in the packages missMDA (for PCA; Husson & Josse, 2016),vegan (for NMDS; Oksanen et al., 2017), lme4 (for LMM;
Bates et al., 2015) and smatr (for OLS; Warton et al., 2012).