Environmental heterogeneity, rather than stability, explains spider assemblage differences between ecosystems

Open ecosystems (e.g., grasslands, prairies, shrublands) tend to be ecologically less stable than closed ones (i.e., forests) and encompass higher spatial heterogeneity in terms of environmental diversity. Such differences are expected to differentially constrain the diversity and structure of the communities that inhabit each of them, but identifying the specific processes driving contrasting biodiversity patterns between open and closed systems is challenging. In order to understand how environmental variability might structure spider assemblages, both between and within open and closed ecosystems, we implement a high throughput multiplex barcode sequencing approach to generate a dataset for 8585 specimens representing 168 species, across open ecosystems within the Canary Islands. Combining these with spider sequences from closed ecosystems within the same islands, we show that spider communities in open ecosystems show higher species richness, higher beta diversity, and higher proportions of rare species but proportionately lower numbers of endemic species than communities in closed ecosystems. We furthermore assess if environmental heterogeneity and habitat stability are the major drivers of such differences by assessing spatial genetic structuring and the influence of bioclimatic variables. Our results point to environmental heterogeneity rather than stability as a major driver of spatial patterns between open and closed ecosystems.

Methods
A total of 25 sites were sampled within open shrublands characterised by the presence of species of Euphorbia across the four western islands of the Canary Archipelago. These data were complemented with an existing dataset including 31 closed laurel forest sites across the same islands. Sampled individuals were examined under a stereomicroscope, and provisional taxonomic assignments were made at the species, genus, or family level. Subsequently, all specimens were DNA-barcoded using the multiplex barcoding approach. Community dissimilarity matrices were estimated at the species level across all sites. Matrices were generated using total β diversity (Sørensen index, βsor) and its additive turnover (Simpson index, βsim) and nestedness (βsne) components. To disentangle the contribution of different factors in the βsor within ecosystems, generalised dissimilarity models (GDM) were applied. A matrix of dissimilarity across sites was implemented as the response variables, and geographical distance and bioclimatic variables from the sampling sites were implemented as predictors. Environmental variables (the same 19 bioclimatic variables defined in WorldClim, Fick, and Hijmans 2017) were downloaded at the finest resolution available (100 m) from CanaryClim (Patiño et al. 2023). Correlations between predictor variables were assessed using Pearson's rank correlation coefficient r, establishing r < |0.7| as a threshold to avoid multicollinearity among the bioclimatic variables used as predictors. Models were then fitted, setting 50 permutations to evaluate the models and variable significance (gdm function, gdm package), to estimate variable importance (gdm.varImp function, gdm package), and to perform deviance partitioning (gdm.partition.deviance function, gdm package). Models were conducted for each ecosystem separately, both across and within islands. Finally, the overall environmental heterogeneity and the spatial structure of this environmental heterogeneity were compared between open and closed ecosystems. First, principal component analyses (PCA) (dudi.pca function, ade4 package; Dray and Dufour 2007) were conducted both at the archipelago and the within-island scale using the 19 bioclimatic variables from Patiño et al. (2023) described before. The betadisper function (vegan package; Oksanen et al. 2022) was used to evaluate the multivariate homogeneity of group dispersions (mean distances to centroid) between ecosystems. To test for temporal stability between ecosystems, generalised linear models (glm function, stats package) were fitted with a Gaussian error structure implementing as dependent variables the following bioclimatic variables related to variation: BIO2 (mean diurnal temperature range), BIO7 (temperature annual range), BIO3 (isothermality), and BIO15 (precipitation seasonality). Finally, multiple regression on distance matrices (MRM function, ecodist package; Goslee and Urban 2007) was implemented to assess the correlation between the geographical and environmental (euclidean) distance matrixes for each ecosystem at both archipelago and island scales. To test for the geographic structuring of genetic variation within PBS, the fixation index NST was calculated using SPAGeDI 1.5. NST was calculated both among islands (each island is considered a population) and within individual islands (each sampling site is considered a population).