A global synthesis reveals biodiversity-mediated benefits for crop production.

Human land use threatens global biodiversity and compromises multiple ecosystem functions critical to food production. Whether crop yield-related ecosystem services can be maintained by a few dominant species or rely on high richness remains unclear. Using a global database from 89 studies (with 1475 locations), we partition the relative importance of species richness, abundance, and dominance for pollination; biological pest control; and final yields in the context of ongoing land-use change. Pollinator and enemy richness directly supported ecosystem services in addition to and independent of abundance and dominance. Up to 50% of the negative effects of landscape simplification on ecosystem services was due to richness losses of service-providing organisms, with negative consequences for crop yields. Maintaining the biodiversity of ecosystem service providers is therefore vital to sustain the flow of key agroecosystem benefits to society.

INTRODUCTION

Natural and modified ecosystems contribute a multitude of functions and services that support human well-being (, ). It has long been recognized that biodiversity plays an important role in the functioning of ecosystems (, ), but the dependence of ecosystem services on biodiversity is under debate. An early synthesis revealed inconsistent results (), whereas subsequent studies suggest that a few dominant species may supply the majority of ecosystem services (, ). It thus remains unclear whether a few dominant or many complementary species are needed to supply ecosystem services.

The interpretation of earlier studies has been controversial because multiple mechanisms underlying changes in ecosystem service response to biodiversity can operate in combination (, ). On one hand, communities with many species are likely to include species responsible for large community-wide effects due to statistical selection. On the other hand, such diverse communities may contain a particular combination of species that complement each other in service provisioning. While these mechanisms imply positive effects of species richness on ecosystem service supply, total organism abundance or dominance of certain species may also drive the number of interactions benefiting ecosystem service supply. Depending on the relative importance of species complementarity, community abundance, and the role of dominant species, different relationships between species richness and ecosystem services can be expected ().

In real-world ecosystems, natural communities consist of a few highly abundant (dominant species) and many rare ones. The importance of richness, abundance, and dominance is likely to be influenced by the extent to which relative abundance changes with species richness () and by differences in the effectiveness and degree of specialization of service-providing communities. However, these three aspects of diversity have typically been tested in isolation and mainly in small-scale experimental settings (, ), while a synthetic study contrasting their relative importance in real-world ecosystems is still lacking. A major limitation to resolving these relationships is a lack of evidence from real-world human-driven biodiversity changes (, ), particularly for ecosystem services in agroecosystems. For instance, changes in richness and total or relative abundance of service-providing organisms in response to land-clearing for agriculture (, ) could alter the flow of benefits to people in different ways compared to experimental random loss of biodiversity.

Over the past half-century, the need to feed a growing world population has led to markedly expanded and intensified agricultural production, transforming many regions into simplified landscapes (). This transformation not only has contributed to enhanced agricultural production but also has led to the degradation of the global environment. The loss of biodiversity can disrupt key intermediate services to agriculture, such as crop pollination () and biological pest control (), which underpin the final provisioning service of crop production (). The recent stagnation or even decline of crop yields with ongoing intensification () indicates that alternative pathways are necessary to maintain future stable and sustainable crop production (–). An improved understanding of global biodiversity-driven ecosystem services in agroecosystems and their cascading effects on crop production is urgently needed to forecast future supplies of ecosystem services and to pursue strategies for sustainable management ().

We compiled an extensive database comprising 89 studies that measured richness and abundance of pollinators, pest natural enemies, and associated ecosystem services at 1475 sampling locations around the world (Fig. 1A). We focused on the ecosystem services of pollination and biological pest control because these services are essential to crop production and have been the focus of much research in recent decades (). We quantified pollinator and pest natural enemy richness as the number of unique taxa sampled from each location (field), abundance as the number of observed individuals, and evenness (or the complementary term, dominance) as the Evar index () that reflects the relative abundance of those taxa. We derived a standardized index of pollination services using measures of pollination success (fruit set or seed set) and plant reproduction (pollen limitation) and of pest control services using measures of natural enemy activity (predation or parasitism rates) and crop damage (pest activity) (see Materials and Methods). We z-transformed each measure separately to remove the effect of differences in measurement scale between indices and inverted values for those measures where low values indicate positive contribution to ecosystem service (i.e., pest activity). We also characterized the 1-km-radius landscape surrounding each field by measuring the percentage of cropland from high-resolution land-use maps. This landscape metric has been used as a relevant proxy for characterizing landscape simplification (, ) and is often correlated with other indicators of landscape complexity (, ).

Fig. 1

Distribution of analyzed studies and effects of richness on ecosystem services provisioning.

(A) Map showing the size (number of crop fields sampled) and location of the 89 studies (further details of studies are given in table S1). (B) Global effect of pollinator richness on pollination (n = 821 fields of 52 studies). (C) Global effect of natural enemy richness on pest control (n = 654 fields of 37 studies). The thick line in each plot represents the median of the posterior distribution of the model. Light gray lines represent 1000 random draws from the posterior. The lines are included to depict uncertainty of the modeled relationship.

Using a Bayesian multilevel modeling approach, we addressed three fundamental yet unresolved questions in the biodiversity-ecosystem function framework: (i) Do richness, abundance, and dominance sustain ecosystem services in real-world ecosystems? (ii) Does landscape simplification indirectly affect ecosystem services mediated by a loss of local community diversity? (iii) How strong are the cascading effects of landscape simplification on final crop production?

RESULTS AND DISCUSSION

We found clear evidence that richness of service-providing organisms positively influenced ecosystem service delivery. This was detected for both pollination (Fig. 1B and table S2) and pest control (Fig. 1C and table S2) and in almost all studies (figs. S1 and S2). As different methods were used in different studies to quantify richness and ecosystem services, we tested the sensitivity of our results to methodological differences. The bivariate relationships between richness and ecosystem services were robust to the taxonomic resolution to which each organism was identified, the sampling methods used to collect pollinators and natural enemies, and the pollination and pest control service measures used (see Materials and Methods).

The positive contribution of richness to service supply was also confirmed in a series of path analyses where we partitioned the relative importance of richness, total abundance, and evenness in driving biodiversity–ecosystem services relationships (Fig. 2). Although richness showed the larger contribution in both models, we also found a direct effect of total abundance (Fig. 2, A and B) and evenness (Fig. 2, C and D) on ecosystem service delivery. While both richness and total abundance showed a positive effect, evenness had a negative effect, suggesting that dominant species contribute more strongly to service supply. This integrative assessment of different aspects of community structure revealed a more multilayered relationship between biodiversity and ecosystem services than has been previously acknowledged. Our results complement previous findings for pollination (, , ) and pest control () and indicate that richness and total and relative abundance are not mutually exclusive but concurrently contribute to support these two key ecosystem services to agriculture. Hence, we find strong support for the role of species-rich communities in supporting pollination and pest control services.

Fig. 2

Direct and indirect effects of richness, total abundance, and evenness on ecosystem services.

(A) Path model of pollinator richness as a predictor of pollination, mediated by pollinator abundance. (B) Path model of natural enemy richness as a predictor of pest control, mediated by natural enemy abundance. (C) Path model of pollinator richness as a predictor of pollination, mediated by pollinator evenness. (D) Path model of natural enemy richness as a predictor of pest control, mediated by natural enemy evenness. Pollination model, n = 821 fields of 52 studies; pest control model, n = 654 fields of 37 studies. Path coefficients are effect sizes estimated from the median of the posterior distribution of the model. Black and red arrows represent positive or negative effects, respectively. Arrow widths are proportional to highest density intervals (HDIs).

Furthermore, we found that landscape simplification indirectly affected ecosystem services by reducing the richness of service-providing organisms. Roughly a third of the negative effects of landscape simplification on pollination were due to a loss in pollinator richness (Fig. 3A and table S5). This effect was even greater for pest control where natural enemy richness mediated about 50% of the total effect of landscape simplification (Fig. 3B and table S5). A similar pattern was also found considering abundance in addition to richness. In this case, landscape simplification indirectly affected ecosystem services by reducing both richness and abundance of service-providing organisms (fig. S4, A and B, and table S6). A consistent richness-mediated effect was also confirmed when we tested the direct and indirect effects of landscape simplification on ecosystem services via changes in both richness and evenness (fig. S4, C and D, and table S7). However, contrary to our expectation, landscape simplification led to higher pollinator evenness. It is likely that landscape simplification has stronger effects on more specialized, rare species and shifts assemblages toward more evenly abundant, mobile, generalist species with a higher ability to use scattered resources in the landscape (, ). These findings, combined with the larger contribution of richness when tested in combination with evenness, suggest that not only dominant species but also relatively rare species contributed to pollination services (, ). The effect of landscape simplification on ecosystem services was minimized when not considering the mediated effect of richness, especially for pest control. We did not find a direct landscape simplification effect on pest control [all highest density intervals (HDIs) overlapped zero; Fig. 3B and table S5]. Together, these results demonstrate strong negative indirect effects of landscape simplification on biodiversity-driven ecosystem services in agroecosystems and the importance of the richness and abundance of service-providing organisms in mediating these effects.

Fig. 3

Direct and indirect effects of landscape simplification on richness of service-providing organisms and associated ecosystem services.

(A) Path model of landscape simplification as a predictor of pollination, mediated by pollinator richness (n = 821 fields of 52 studies). (B) Path model of landscape simplification as a predictor of pest control, mediated by natural enemy richness (n = 654 fields of 37 studies). Path coefficients are effect sizes estimated from the median of the posterior distribution of the model. Black and red arrows represent positive and negative effects, respectively. Arrow widths are proportional to HDIs. Gray arrows represent nonsignificant effects (HDIs overlapped zero).

Last, for a subset of the data that had crop production information (674 fields of 42 studies), we found that the cascading effects of landscape simplification, mediated through richness and associated ecosystem services, led to lower crop production. This was detected for both pollination (Fig. 4A and table S8) and pest control (Fig. 4B and table S8). Specifically, landscape simplification reduced both pollinator and natural enemy richness, which had indirect consequences for pollination and pest control and, in turn, decreased crop production. For pest control, a positive link with crop production was detected in fields where the study area was not sprayed with insecticides during the course of the experiment (Fig. 4B) but not when considering all sites combined (with and without insecticide use; fig. S4). In sprayed areas, we did not find a pest control effect (all HDIs overlapped zero), probably because effects were masked by insecticide use (, ), indicating that insecticide use undermines the full potential of natural pest control. A positive link with crop production was detected, although measures used to estimate pest control (natural enemy and pest activity) were not direct components of crop production, as was the case of pollination measures (fruit or seed set). We also found an indirect influence of dominance (Fig. 4 and table S8) but not of abundance (fig. S4 and table S9) on crop production. However, the indirect dominance pathway was weaker than that via richness. An abundance effect may not have been detectable because of the lower sample size in these submodels with crop production. Although only available from a subset of the data, this result supports the hypothesis that the effects of landscape simplification can cascade up to reducing the final provisioning service of crop production.

Fig. 4

Direct and cascading effects of landscape simplification on final crop production via changes in richness, evenness, and ecosystem services.

(A) Path model representing direct and indirect effects of landscape simplification on final crop production through changes in pollinator richness, evenness, and pollination (n = 438 fields of 27 studies). (B) Path model representing direct and indirect effects of landscape simplification on final crop production through changes in natural enemy richness, evenness, and pest control [only insecticide-free areas were considered in the model (n = 185 fields of 14 studies)]. Path coefficients are effect sizes estimated from the median of the posterior distribution of the model. Black and red arrows represent positive and negative effects, respectively. Arrow widths are proportional to HDIs. Gray arrows represent nonsignificant effects (HDIs overlapped zero).

Our findings suggest that some previously inconsistent responses of natural enemy abundance and activity to surrounding landscape composition () can be reconciled by considering richness in addition to abundance. Although richness and abundance are often correlated, their response to environmental variation can differ. This was evident in the path analysis showing a strong effect of landscape simplification on richness but only a marginal effect on abundance (fig. S3B). Moreover, effect sizes for natural enemy abundances in individual studies (fig. S5) showed similarly inconsistent responses to the previous synthesis (). Likewise, results for natural enemy activity showed inconsistent responses to surrounding landscape composition in both syntheses [Fig. 3B and ()].

Using an integrative model to assess key ecological theory, we demonstrate that the negative effects of landscape simplification on service supply and final crop production are primarily mediated by loss of species. We found strong evidence for positive biodiversity–ecosystem service relationships, highlighting that managing landscapes to enhance the richness of service-providing organisms (, ) is a promising pathway toward a more sustainable food production globally. In an era of rapid environmental changes, preserving biodiversity-driven services will consistently confer greater resilience to agroecosystems, such that we could expect improved crop production under a broader range of potential future conditions.

MATERIALS AND METHODS

Database compilation

We compiled data from crop studies where measures of richness and abundance of service-providing organisms (pollinators or natural enemies) and associated ecosystem services (pollination and biological pest control) were available for the same sites. If available, we also included information on yield. Studies were identified by first searching the reference lists of recent meta-analyses (, , , , ) and then directly contacting researchers. For pest control, data were mostly provided from a recent pest control database (). Of 191 researchers initially contacted, 86 provided data that met our criteria. As similar studies were frequently performed in the same area, occasionally in the same year, and studies with multiple years usually used different sites each year, we did not nest year within study. Likewise, some studies collected data in different crops. Thus, we considered each year of multiyear studies (that is, 10 studies) and each crop of multicrop studies (that is, 5 studies) to be an independent dataset and used study-crop-year combinations as the highest hierarchical unit. Overall, we analyzed data from 89 studies and 1475 fields in 27 countries around the world (table S1). Twenty-nine crops were considered, including a wide array of annual and perennial fruit, seed, nut, stimulant, pulse, cereal, and oilseed crops. Studies represented the spectrum of management practices, including conventional, low-input conventional, integrated pest management, and organic farming. Few studies compared different management practices in similar landscapes, crop types, or sampled year. In 76% of fields, pest control experiments were performed in insecticide-free areas. In some fields, this information was not available (7%) or insecticides were applied (17%).

Pollinator and pest natural enemy richness and abundance

Studies used a broad range of methods, which we categorized as active or passive (), to sample pollinators or natural enemies. Active sampling methods included netting pollinators seen on crop flowers, hand-collecting individuals on plants, observational counting, sweep-netting, and vacuum sampling. Passive sampling methods were malaise traps, pan traps, pitfall traps, and sticky cards. Active sampling was performed in 85% of pollinator sampling fields and in 50% of natural enemy sampling fields.

Pollinators included representatives from the orders Hymenoptera, Diptera, Lepidoptera, and Coleoptera. Bees (Hymenoptera: Apoidea) were the most commonly observed pollinators and included Apis bees (Apidae: Apis mellifera, Apis cerana, Apis dorsata, and Apis florea), stingless bees (Apidae: Meliponini), bumblebees (Apidae: Bombus spp.), carpenter bees (Apidae: Xylocopini), small carpenter bees (Apidae: Ceratinini), sweat bees (Halictidae), long-horned bees (Apidae: Eucerini), plasterer bees (Colletidae), mining bees (Andrenidae), and mason bees (Megachilidae). Non-bee taxa included syrphid flies (Diptera: Syrphidae), other flies (Diptera: Calliphoridae, Tachinidae, and Muscidae), butterflies and moths (Lepidoptera), various beetle families (Coleoptera), and hymenopterans including ants (Formicidae) and the paraphyletic group of non-bee aculeate wasps. Natural enemies included ground beetles (Coleoptera), flies (Diptera), spiders (Aranea), hymenopterans including ants (Formicidae) and wasps, bugs (Hemiptera), thrips (Thysanoptera), net-winged insects (Neuroptera), bats, and birds.

We calculated pollinator and natural enemy richness as the number of unique taxa sampled per study, method, and field. A taxon was defined as a single biological type (that is, species, morphospecies, genus, and family) determined at the finest taxonomic resolution to which each organism was identified. In almost 70% of cases, taxonomic resolution was to species level (averaged proportion among all studies), but sometimes, it was based on morphospecies (15%), genus (8%), or family levels (7%). Taxon richness per field varied between 1 and 49 for pollinators and between 1 and 40 for natural enemies. Abundance reflected the sum of individuals sampled per study, method, and field. Pollinator richness was calculated either including or excluding honey bees (A. mellifera). A. mellifera was considered as the only species within the honey bee group for consistency across all datasets (). Other Apis bees (that is, A. cerana, A. dorsata, and A. florea) were not pooled into the honey bee category as the large majority of observed individuals were derived from feral populations. Feral and managed honey bees were analyzed as a single group because they cannot be distinguished during field observations. Feral honey bees were uncommon in most studies except for those in Africa and South America (A. mellifera is native in Africa, while it was introduced to the Americas). In studies with subsamples within a field (that is, plots within fields or multiple sampling rounds within fields), we calculated the total number of individuals and unique taxa across these subsamples.

We also calculated evenness using the Evar index (). This index is based on the variance of log abundances (centered on the mean of log abundances) and then appropriately scaled to cover 0 to 1 (0, maximally uneven and 1, perfectly even).

Pollination and pest control services

As different methods were used to quantify pollination or pest control services across studies, standardization was necessary to put all the indices on equivalent terms. Therefore, we transformed each index y in each field i in each study j using z scores. We preferred the use of z scores over other transformations (for example, division by the maximum) because z scores do not constrain the variability found in the raw data, as do other indices that are bounded between 0 and 1. We used the proportion of flowers that set fruit (that is, fruit set), the average number of seeds per fruits (that is, seed set), or the estimated measures of pollinator contribution to plant reproduction (that is, differences in fruit weight between plants with and without insect pollination, hereafter Δ fruit weight) as measures of pollination services. We then converted these measures into the pollination index. The pest control index was calculated using measures of natural enemy activity or pest activity. Natural enemy activity was measured by sentinel pest experiments where pests were placed in crop fields and predation or parasitism rates were monitored or by field exclosure experiments where cages were used to exclude natural enemies to quantify differences in pest abundance or crop damage between plants with and without natural enemies. Pest activity was measured as the fraction or amount of each crop consumed, infested, or damaged. We inverted standardized values of pest activity by multiplying by −1, as low values indicate positive contributions to the ecosystem service.

Crop production

Depending on the crop type, marketable crop yield is valued by farmers not only in terms of area-based yield but also in terms of fruit or seed weight [for example, in coffee, sunflower, or strawberry fields; (, )] or seed production per plant [for example, in seed production fields; ()]. Moreover, area-based yield and within-plant yield are often correlated (, ). Thus, we used both area-based yield and within-plant yield as measures of final crop production. Within-plant yield was measured by the total number (or mass) of seeds or fruits per plant or by fruit or seed weight. In addition, in this case, we standardized variables (z scores) to put all the indices on equivalent terms.

Landscape simplification

Landscapes were characterized by calculating the percentage of cropland (annual and perennial) within a radius of 1 km around the center of each crop field. This landscape metric has been used as a relevant proxy for characterizing landscape simplification (, ) and is often correlated with other indicators of landscape complexity (, ). Moreover, we used this metric because cropland data are readily accessible from publicly available land cover data and are more accurate than other land use types such as forests and grasslands (), especially when detailed maps are not available. The 1-km spatial extent was chosen to reflect the typical flight and foraging distances of many insects including pollinators (, ) and natural enemies (, ). For studies where this information was not supplied by the authors, land uses were digitized using GlobeLand30 (), a high-resolution map of Earth’s land cover. The derived land cover maps were verified and, if necessary, corrected using a visual inspection of satellite images (Google Earth). We then calculated the percentage of cropland within the radius using Quantum GIS 2.18 (Open Source Geospatial Foundation Project; http://qgis.osgeo.org). The average percentage of cropland was 67.5% for pollination studies and 41.5% for natural enemy studies.

Data analysis

Data standardization

Before performing the analyses, we standardized the predictors (abundance, richness, and landscape simplification) using z scores within each study. This standardization was necessary to allow comparisons between studies with differences in methodology and landscape ranges (). Moreover, this allows the separation of within-study effects from between-study effects. This separation is important because it prevents the risk of misinterpreting the results based on studies differing in methodology and landscape gradients, by erroneously extrapolating from between-study effects to within-study effects or vice versa (, ).

Relationship between richness and ecosystem services

The relationship between richness of service-providing organisms and related ecosystem services (Fig. 1, B and C) was estimated from a Bayesian multilevel (partial pooling) model that allowed the intercept and the slope to vary among studies (also commonly referred to as random intercepts and slopes), following the equationwhere ES is the ecosystem service index (pollination or pest control depending on the model), RIC is richness of service-providing organisms (pollinator or natural enemy richness depending on the model), and j[ represents observation i of study j. This partial-pooling model estimated both study-level responses [yielding an estimate for each study (β)] and the distribution from which the study-level estimates were drawn, yielding a higher-level estimate of the overall response across crop systems (μβ). In addition, it accounted for variation in variance and sample size across observations (for example, studies). The intercepts α and slopes β varied between studies according to a normal distribution with mean μ and SD σ. Independent within-study errors also followed a normal distribution ε ~ N(0,σ). We used weakly informative priors: normal (0, 10) for the population-level parameters (α, β) and half–Student t (3, 0, 5) for the group-level SD and residual SD.

Direct and indirect effects of richness, abundance, and evenness on service provisioning

As natural communities vary not only in number of species but also in relative abundance of each species (evenness) and the total number of individuals (abundance), it is important to incorporate these attributes when assessing or modeling biodiversity effects (, ). In a Bayesian multivariate response model, a form of path analysis, we partitioned the relative importance of richness and total and relative abundance in driving biodiversity-ecosystem relationships. We hypothesized that richness drives both abundance and evenness according to a revised version of the “more individuals hypothesis” () that proposes that more species can exploit more diverse resources and may therefore maintain more individuals than species-poor communities. Specifically, we tested (i) whether richness per se directly influences ecosystem services or is instead mediated by abundance (Eq. 1, 3; Fig. 2, A and B) and (ii) whether richness per se directly influences ecosystem services or is instead mediated by evenness (Eq. 2, 4; Fig. 2, C and D). Before analysis, we also checked for data collinearity among abundance, evenness, and richness by calculating the variance inflation factor (VIF). No signal of collinearity was detected in either model (VIFs were below 1.8). In a preliminary analysis, we also tested a possible interaction between richness and evenness. No significant interaction was found in both pollination and pest control models. To illustrate, we first show the univariate multilevel (partial pooling) models following these equationswhere RIC is richness of service-providing organisms, ABU is abundance, EVE is evenness, ES is the ecosystem service index, and the index j[ represents observation i of study j. We then specified both multivariate multilevel models in a matrix-vector notion (), as followswhere Y is the matrix of response variables with observations i as rows and variables r as columns, X is the matrix of all predictors for response r, B are the regression parameters (α and β) for response r, M represents the mean of the distribution of the regression parameters, and ∑ is the covariance matrix representing the variation of the regression parameters in the population groups. We used weakly informative priors: normal (0, 10) for the population-level parameters (α, β) and half–Student t (3, 0, 5) for the group-level SD and residual SD. In building the model, we ensured that no residual correlation between ES and RIC, between ES and EVE, or between ES and ABU was estimated [see the “set_rescor” function in the package brms; ()].

Direct and indirect effects of landscape simplification on ecosystem services

To estimate the direct and indirect effects of landscape simplification on richness and associated ecosystem services, we used two models. In a Bayesian multivariate response model with causal mediation effects (hereafter, mediation model), we tested whether landscape simplification directly influences ecosystem services or is mediated by richness. Mediation analysis is a statistical procedure to test whether the effect of an independent variable X on a dependent variable Y (X → Y) is at least partly explained via the inclusion of a third hypothetical variable, the mediator variable M (X → M → Y) (). The three causal paths a, b, and c′ correspond to X’s effect on M, M’s effect on Y, and X’s effect on Y accounting for M, respectively. The three causal paths correspond to parameters from two regression models, one in which M is the outcome and X is the predictor and one in which Y is the outcome and X and M are the simultaneous predictors (fig. S6). From these parameters, we can compute the mediation effect (the product ab, also known as the indirect effect) and the total effect of X on Y

Thus, the total causal effect of X, which is captured by the parameter c, can be decomposed precisely into two components, a direct effect c′ and an indirect (mediation) effect ab (the product of paths a and b). The model included the ecosystem service index as response, landscape simplification as predictor, and richness as mediator (Fig. 3). The separate regression models that made up the Bayesian multivariate multilevel model followed these equations

We then compiled a multilevel path analysis testing the direct and indirect effects of landscape simplification on ecosystem services via changes in both richness and abundance (Eq. 5, 6, 8; fig. S4, A and B) or richness and evenness (Eq. 5, 7, 9; fig. S4, C and D). The separate regression models that made up the model followed these equationswhere RIC is richness of service-providing organisms, LAND is landscape simplification measured as the percentage of arable land surrounding each study site, ABU is abundance, EVE is evenness, ES is the ecosystem service index, and the index j[ represents observation i of study j. We then specified multivariate multilevel models in a matrix-vector notion, as explained above. The mediation analysis was implemented using the R package sjstats [v. 0.15.0; ()].

Cascading effects of landscape simplification on final crop production

For 42 studies and 675 fields (pollination model, n = 438 fields of 27 studies; pest control model, n = 236 fields of 15 studies; table S1), the data allowed us to use a multilevel path analysis to examine cascading effects of landscape simplification on final crop production via changes in richness, abundance, evenness, and ecosystem services. We expected that (i) landscape simplification would have a direct effect on richness, abundance, and dominance of service-providing organisms and (ii) richness and abundance dominance would relate positively to intermediate services, which, in turn, (iii) would increase final crop production (Fig. 4 and fig. S4). The indirect effects of richness and abundance (Eq. 10, 11, 13, 15; Fig. 4) or richness and evenness (Eq. 10, 12, 14, 15; fig. S4) were tested separatelywhere RIC is richness of service-providing organisms, LAND is landscape simplification measured as the percentage of arable land surrounding each study site, ABU is abundance, EVE is evenness, ES is the ecosystem service index, PROD is crop production, and the index j[ represents observation i of study j. We specified a multivariate multilevel model in a matrix-vector notion, as explained above.

Parameter estimation

All analyses were conducted in the programming language Stan through R (v. 3.4.3) using the package brms [v. 2.2.0; ()]. Stan implements Hamiltonian Monte Carlo and its extension, the No-U-Turn Sampler. These algorithms converge much more quickly than other Markov chain Monte Carlo algorithms especially for high-dimensional models (). Each model was run with four independent Markov chains of 5000 iterations, discarding the first 2500 iterations per chain as warm-up and resulting in 10,000 posterior samples overall. Convergence of the four chains and sufficient sampling of posterior distributions were confirmed by (i) the visual inspection of parameter traces, (ii) ensuring a scale reduction factor below 1.01, and (iii) an effective size (neff) of at least 10% of the number of iterations. For each model, posterior samples were summarized on the basis of the Bayesian point estimate (median), SE (median absolute deviation), and posterior uncertainty intervals by HDIs, a type of credible interval that contains the required mass such that all points within the interval have a higher probability density than points outside the interval (). The advantage of the Bayesian approach is the possibility to estimate not only expected values for each parameter but also the uncertainty associated with these estimates (). Thus, we calculated 80, 90, and 95% HDIs for parameter estimates.

Sensitivity analyses

Given that different methods were used in different studies to quantify richness, ecosystem services, and final crop production, we measured the sensitivity of our results to methodological differences.

(1) We verified whether treating each annual dataset from multiyear studies separately could incorrectly account for the dependence of the data. We refitted the model testing the relationship between richness and ecosystem services including year nested within study (that is, study defined as study-crop combination). Then, we compared models (year-independent model versus year-nested model) using leave-one-out (LOO) cross-validation, a fully Bayesian model selection procedure for estimating pointwise out-of-sample prediction accuracy (). We calculated the expected log pointwise predictive density , using the log-likelihood evaluated at the posterior simulations of the parameter values. Model comparison was implemented using R package loo [v. 2.0.0; ()]. We found that the year-nested model had a lower average predictive accuracy than the year-independent model for both pollination and pest control and therefore retained the year-independent model in our analysis.

(2) We verified whether taxonomic resolution influenced the interpretation of results. We recalculated richness considering only organisms classified at the fine taxonomy level (species or morphospecies levels) and refitted the model testing the effect of richness on ecosystem services. We found no evidence that taxonomic resolution influenced our results. With a fine taxonomic resolution, the effects of richness on ecosystem services (βpollinators = 0.1535, 90% HDIs = 0.0967 to 0.2141; βenemies = 0.2264, 90% HDIs = 0.1475 to 0.3065; table S2) were nearly identical to the estimates presented in the main text (βpollinators = 0.1532, 90% HDIs = 0.0892 to 0.2058; βenemies = 0.2093, 90% HDIs = 0.1451 to 0.2779; table S2).

(3) We verified whether the sampling methods used to collect pollinators (active versus passive sampling techniques) influenced the relationship between pollinator richness and pollination using Bayesian hypothesis testing (). Passive methods do not directly capture flower visitors and may introduce some bias (for example, they may underestimate flower visitors). However, our estimate was not influenced by sampling method (the one-sided 90% credibility interval overlapped zero; table S11). In accordance with this finding, the evidence ratio showed that the hypothesis tested (that is, estimates of studies with active sampling > estimates of studies with passive sampling) was only 0.78 times more likely than the alternative hypothesis.

(4) We verified whether methodological differences in measuring pollination and pest control services influenced the relationship between richness and ecosystem services. Using Bayesian hypothesis testing, we tested whether the estimates differed among methods. The two-sided 95% credibility interval overlapped zero in all comparisons (estimates did not differ significantly; table S12), indicating that our estimate was not influenced by methodological differences in measuring ecosystem services. Furthermore, we tested effects including only inverted pest activity as a reflection of pest control. We found positive effects of natural enemy richness on inverted pest activity (β = 0.1307, 90% HDIs = 0.0102 to 0.2456), indicating that results were robust to the type of pest control measure considered.

(5) As honey bees are the most important and abundant flower visitors in some locations, we verified the potential influence of honey bees on our results by refitting the path models testing direct and indirect effects of richness, abundance, and evenness on pollination services with honey bees. A positive direct contribution of richness to pollination was confirmed even after including honey bees (fig. S7). In these models, both abundance and evenness showed a larger effect compared to models without honey bees (Fig. 2).

(6) Insecticide application during the course of the experiment could mask the effect of pest control on crop production (, ). We verified the potential influence of insecticide application on our results by refitting the model considering only fields where the study area was not sprayed with insecticide during the course of the experiment (n = 85 fields of 14 studies). We found a pest control effect that was masked when considering all sites combined (with and without insecticide; fig. S4). We therefore show the insecticide-free model in the main text (Fig. 4).

(7) We verified the consistency of our results considering only studies that measured area-based yield (submodel). Only the cascading effects of landscape simplification on final crop production via changes in richness were tested in a simplified model. We found no evident differences between the submodel (fig. S8) and the full model presented in the main text (Fig. 4).