Terrestrial land-cover type richness is positively linked to landscape-level functioning

Biodiversity–ecosystem functioning (BEF) experiments have shown that local species richness promotes ecosystem functioning and stability. Whether this also applies under real-world conditions is still debated. Here, we focus on larger scales of space, time and ecological organization. We develop a quasi-experimental design in which we relate land-cover type richness as measure of landscape richness to 17-year time series of satellite-sensed functioning in 4974 landscape plots 6.25 or 25 ha in size. We choose plots so that landscape richness is orthogonal to land cover-type composition and environmental conditions across climatic gradients. Landscape-scale productivity and temporal stability increase with landscape richness, irrespective of landscape plot size. Peak season near-infrared surface albedo, which is relevant for surface energy budgets, is higher in mixed than in single land-cover type landscapes. Effect sizes are as large as those reported from BEF-experiments, suggesting that landscape richness promotes landscape functioning at spatial scales relevant for management.


Pairwise interactions among land-cover units
We analyzed pairwise interactions of land-cover units by mechanistic diallel analysis 1 . measures the deviation from this expectation, i.e. the net interaction of a and b. The diallel framework thus allows to identify positive or negative interactions between two components of a unit, similar to the additive partitioning method of Loreau & Hector 2 . An advantage of the diallel method is that it can be applied to cases where data is available at the whole unit level but not separately for the individual components a and b. A further difference to the additive partitioning is that SCAs are corrected for differences in the average performance of each unit across combinations, which is not the case for the NE and CE in the additive partitioning method.

4
We determined GCA and SCA values by fitting a mechanistic diallel model to the subset of our data where landscape richness ≤ 2, separately for every block.
Conceptually, single land-cover landscape plots were treated as combination of twice the same land-cover type. The model matrix was constructed by superimposing effects of the two component land-cover types (i.e. both GCA terms) using the and() special function of ASReml (VSN International, Hemel Hemsted, UK; see ref. 3 ).
For the subsequent analyses, we used these block-wise GCA and SCA values as new dependent variables. We tested whether SCAs differed between single land-cover landscapes and mixed landscapes using a model with landscape richness as fixed effect.
We then tested whether GCAs depended on land-cover type. Similarly, we tested whether SCAs depended on land-cover pairs and, if significant, performed pair-wise multiple comparisons tests (Tukey's HSD, function glht in R-library multcomp 4 ).

Correlation of landscape richness with species richness
We tested for correlations between local species richness and land-cover type richness (i.e. landscape richness) in our study area using data from the Swiss Biodiversity Monitoring Program (BDM; biodiversitymonitoring.ch; refs. 5,6 ). The presence of vascular plant species within ecosystems ("Z9 Indicator"; ref. 6 ) and within landscapes ("Z7 Indicator"; ref. 6 ) was recorded in five-year intervals since the year 2001. The Z9 data were recorded in over 1400 circular 10m 2 plots. The Z7 data was recorded along two transects across more than 500 landscape plots of 1 km 2 in size. Z7 and Z9 plots both are regularly spread across Switzerland.

5
We calculated herbaceous α-species richness in the Z9 plots classified as meadow or another form of grassland (Sgrass) for the 2001-2016 period. The 10-m 2 Z9 plots were too small to reasonably assess the species richness of woody species, including trees, and we therefore used woody species richness along the Z7 transects (Sforest). In principle, such a transect could cross different land-cover types, but woody species occur primarily in forest and the land-cover type forest varies mostly at larger environmental and spatial scales. We therefore consider these data a reasonable approximation of α-richness of woody forest species.
We then determined landscape richness for the Z7 and Z9 plots. For the 10-m 2 Z9 plots, we calculated landscape richness in quadratic areas 250×250 m and 500×500 m in size with the Z9 plot in their center. For the 1 km 2 Z7 plots, we calculated landscape richness for a 250×250 m (and 500×500 m) area basis by dividing the Z7 plots into 16 (and 4) quadrats and averaging the landscape richness values obtained for these subquadrats.
Next, we removed all Z9 plots and Z7 sub-quadrats that we would not have selected in our main study design because (1) the evenness of land-cover types did not comply with the selection criterion used in the main study, (2) the plot contained compositions that were not included in the main data set (e.g. the land-cover type "urban green"), (3) plots had landscape richness >4, (4) plots had altitude, slope inclination, or north aspect of the slope values outside the range used in the main analysis, or (5) plots were in blocks that were not part of the main analysis (see Figure 1b in main text). Partial correlation coefficients between landscape richness and Z9 or Z7 species richness were calculated after adjusting for block effects (Supplementary Table 1).

Long-range interactions among land-cover types
Functioning metrics were only available at the level of entire landscape plots.
Complementing the diallel analysis, we tested for effects of surrounding land-cover types on landscape plots containing a single land-cover type (landscape richness=1). Forests were the only such landscapes that were dominated by vegetation and occurred in enough replicates (n=390 and 131 for 250×250 m and 500×500 m landscape plots, respectively). We characterized the surrounding of these forest landscapes by the fractional amount of each of the eight aggregated land-cover types (see Methods in the main text). For 250×250 m landscape plots, we focused on a circular area with radius 500 m from the center of the landscape plot. For 500×500 m landscape plots, this radius was increased to 1000 m. We tested for effects of the fractional cover of each surrounding land-cover type using a linear model with block followed by the fractional cover of the land-cover type under consideration.

Spatial configurational diversity of land-cover units
We determined metrics characterizing the spatial configurational diversity of land-cover units (Supplementary Table 2) within all 250×250 m and 500×500 m landscape plots.
Next, we determined unadjusted and block-adjusted partial Pearson correlations between these variables, including our design variable log(LR) (i.e. log-transformed landscape richness, see Methods, main text) and the landscape functioning variables 7 productivity, albedo and their temporal stability. We then tested the effects of configurational diversity on the landscape functioning variables by using these metrics instead of log(LR) (models described in Methods, main text). Finally, we tested whether configurational diversity metrics could explain variation in landscape functioning variables if fitted after log(LR) in the statistical models, i.e. whether configurational diversity could explain additional variation not already explained by log(LR).

Pairwise interactions among land-cover units
Neither GCAs nor SCAs differed between 250×250 m and 500×500 m landscape plots In accordance with our analysis of the main text, SCAs were significantly higher in mixed compared to single land-cover type landscape plots for the primary productivity variables and their temporal stability (EVI ̅̅̅̅̅ , EVIGS,CV EVI ̅̅̅̅̅ −1 , and CV EVI GS −1 ) for the 250×250 m landscapes (all F1,19≥5.8; P≤0.026). These effects were similar but less significant in the 500×500 m landscapes.

Correlation of landscape richness with species richness
For both 250×250 m and 500×500 m landscapes we found that Person correlation coefficients between landscape richness and species richness of grassland and forest plants (Sgrass and Sforest, respectively) were relatively small (Supplementary Table 1).

Long-range interactions among land-cover types
Testing

Spatial configurational diversity of land-cover units
The different configurational diversity metrics were highly correlated across and within