Biodiversity-ecosystem functioning (BEF) experiments have established generally positive species richness-productivity relationships in plots of single ecosystem types. Here, we analyzed effects of landscape-level diversity, measured as the number of land-cover types (different ecosystems) per 250 × 250 m, across all of North America. We find that this metric is positively related to landscape-wide remotely-sensed primary production, and that a higher number of land-cover types also is associated with greater temporal stability of productivity, and with accelerated 20-year greening trends, in particular at high latitudes. Species diversity was correlated with landscape-level productivity, but the effect of species diversity and landscape diversity were independent. This indicates that diversity-functioning patterns resembling the ones at smaller scales also exist at higher levels of biological organization.

Data was collected by processing satellite-remote sensing products collected with the MODIS instrument, at 250m pixel resolution. Land-cover type information was extracted at 30-m spatial resolution from the Commission for Environmental Cooperation’s North American Land Monitoring System’s map (CEC map, based on Landsat-7 satellite imagery), and from the global GlobeLand30 map (GLC map, based on Landsat-5 and China Environmental Disaster Alleviation Satellite (HJ-1) imagery). We focused on the land covers forest, grassland, shrubland, agriculture, wetland and urban, combining the different forest types distinguished in the CEC map.

Study plots were selected across North America to form a quasi-experimental study design with 3x6° latutude x longitude tiles that were further divided into 16 ecoregions. Within each tile x ecoregion combination, parallel experimental sub-designs spanning gradients in land-cover type richness were formed. Plots were selected so that land-cover type richness was orthogonal with the average area fraction of each land-cover type found at each richness level, and so that the richness gradient was orthogonal with important environmental factors such as altitude.

As productivity metric, we used MODIS EVI indices from the Terra satellite (years 2000-2019) and fitted harmonic time series to the data based on Fourier synthesis to model annual phenology curves.