Drivers of species taxonomic richness across beaches and levels
We tested whether number of species was significantly different across beach levels using generalized least squares models (GLS), which allow us to account for the potential effect of various spatial autocorrelation structures, using the ‘gls’ function in the package nlme version 3.1-153 (Piñeiro et al., 2007; Supplementary methods 1.4). In each model, we included the number of species at each sample as the response variable and beach level, mean grain size, beach length, as well as all bioclimatic and oceanographic variables as predictors, in addition to a spatial structure. Model assumptions and model fit were assessed for each model with the R package performance v. 0.7.3 (Lüdecke et al., 2021), testing for distribution of residuals, homoscedasticity, multicollinearity, and influential observations (Zuur et al. 2010). For the models that included a set of predictors with both categorical and continuous variables, which are not easily summarised as output of GLS, we used Type II ANOVA tables using the function ‘Anova’ in the R package car v. 3.0.10 (Fox and Weisberg, 2019).