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