Vegetation and vertebrate abundance as drivers of bioturbation patterns along a climate gradient

classified as any plant type (herbs, shrubs, cacti, trees) to the amount of all pixels. Following [42], the average elevation (hillside elevation) and the hillslope of each plot were estimated based on high resolution Lidar data [43]. Statistical analyses For the burrows of invertebrates and vertebrates we analyzed the allometric relationship between their depth (de) and their diameter (mean of dv and dh). We regressed the mean diam- eter of the entrance versus the depth using the log10-transformed values of both variables and then determining the slope. In an isometric relationship, the log-transformed variables should be linearly related to a slope of one [44]. Since diameter and depth were measured with roughly equal error, in addition to an ordinary least squares (OLS) regression, we estimated the slope using a reduced major axis (RMA) regression [45]. To assess a deviation from a slope of one, we used the offset argument available in most regression functions. With the diameter serving as the independent variable and the offset, the estimate tests for deviations from one. For the slope of the RMA regression, we used the standard error and a t-test to test for devia- tions from one. The same approach was applied to the regression between excavated soil vol- ume and hole density. To analyze the interaction of abiotic and biotic components in bioturbation activity and quantity, we applied generalized linear mixed effect models (GLMMs). We used hole density or excavated soil volume as response variables, site, season, hillside elevation and hillslope as

abiotic fixed predictors and vegetation cover and animal group as biotic fixed predictors. The study plots were used as a random factor (Table 1). All data of the GLMM parameters can be obtained from S7 Table in S4 Appendix. We also included interaction terms between site and all other fixed predictor variables and between vegetation cover and taxon. We standardized the fixed predictor hillside elevation for each site because it varied and could not be assigned separately to each of the sites. We performed GLMMs for the 12 plots within each site (total of 48 plots) in the first field campaign, conducted in the southern-hemispheric autumn, and in the 20 plots within each site (total of 80 plots) during the second campaign conducted in the southern-hemispheric spring. Separation of the hole density of invertebrates and vertebrates resulted in 256 measurements (2 × (48 + 80)). For the GLMM of the excavated soil volume, we log10-transformed data for hole density and excavated soil volume to achieve normality of the residuals. For the log10-transformation, we only considered plots with a hole density > 0 [no/ 100 m2]. Thus, 46 plots without holes were not included in the GLMM for excavated soil vol- ume, resulting in 210 valid measurements. Additionally, we integrated the interaction between hole density and taxon as another fixed predictor. All statistical analyses were performed using the R statistical environment (version 1.3.1093). We used the lmodel2 package [46] for OLS and RMA regression analysis. For the

GLMM, we employed the buildmer function [47] of the lme4 package [48] to perform back- ward stepwise selection. To determine the proportion of variation explained by the model in Table 1. Summary of all variables used in the GLMM. Response variable Abiotic fixed predictors Biotic fixed predictors Random factor Hole density or excavated soil volume Site Vegetation cover Plot number Season Hillside elevation Animal group Hillslope Depicted are the response variables, fixed predictors (abiotic and biotic) and the random factor. https://doi.org/10.1371/journal.pone.0264408.t001 PLOS ONE Vegetation and vertebrate abundance as drivers of bioturbation patterns along a climate gradient PLOS ONE | https://doi.org/10.1371/journal.pone.0264408 March 4, 2022 4 / 14