Quantifying crop pollinator dependence and its heterogeneity using multi-level meta-analysis

1. Biotic pollination can benefit crop production, but its effects are highly variable. To maximise benefits this ecosystem service, we the factors variation that ecological intensification We focus on understanding the benefits of pollination to faba bean Vicia faba . We use a literature review followed by multi-level meta-analysis to estimate overall benefits of pollination to faba bean yield and to quantify variation (heterogeneity) these benefits associated with different contextual factors (e.g. genotype, growing


| INTRODUC TI ON
The contribution of biotic pollination to crop yield has been recently valued at $235-577 billion globally each year .
Understanding this ecosystem service can help us to maximise those pollination benefits by guiding effective decision-making in policy, crop breeding and agronomy. This ecological intensification offers a more sustainable route to increased crop productivity than the use of conventional inputs and agrochemicals .
However, there is significant variation in the contribution of pollinators to crop yield; it varies around the world due to differences in the crop types cultivated , the diversity and abundance of pollinator communities (Woodcock et al., 2019), and many other contextual factors (Garibaldi et al., 2018;Tamburini et al., 2019). It is essential to quantify and explain this variation so that ecological intensification can be more optimally applied in specific contexts. With a better understanding, pollination could increasingly be considered an agronomic input that can complement or replace conventional inputs (Garratt et al., 2018;Tamburini et al., 2017) with more predictable returns for investments in on-farm ecological management . We advocate and demonstrate the use of multilevel meta-analysis to understand variability in pollination benefit and predict the range of pollination dependence across different contexts.
Isolated experimental and observational studies are inadequate to assess and understand the benefits of biotic pollination to a crop.
The primary benefit of biotic pollination to plants is via increased ovule fertilisation, and an increased number of seeds for a plant to develop. If a plant has sufficient resources, there can be a corresponding increase in crop yield (Garibaldi et al., 2018). High soil nutrient availability can therefore make plants more responsive to pollination (Garratt et al., 2018). Conversely, adverse weather or herbivory can suppress the benefits of pollination by damaging floral organs and developing seeds, or by causing resource limitation Lundin et al., 2013). Each isolated study represents a specific combination of nutrient availability, weather and pest pressure, all of which limit the applicability of the findings to different contexts (Garibaldi et al., 2018;Tamburini et al., 2019). Quantitative research synthesis, using 'review-generated' evidence (Cooper, 2015), is therefore needed to understand and account for the variation associated with these different contextual factors.
To encourage ecological intensification by farmers and reap its associated environmental benefits, it is important to report crop responses in economically relevant terms . The use of different metrics to measure crop yield responses to pollination across the literature therefore presents a challenge. Some studies may be concerned only with fertilisation and therefore present abstract fertilisation metrics such as the number of seeds or pods per flower (Drayner, 1959). Due to time or labour constraints, other studies may forgo yield mass measurement and only record simpler metrics such as the number of seeds in a representative sample of pods. The relationships between different yield components can be nonlinear  and as discussed above, contextual factors can modify the translation of pollination into crop yield and economic output. A better understanding of the relationships between different yield metrics could help practitioners to generalise across findings and apply literature from more relevant contexts.
To determine the effect of biotic pollination, it is typically necessary to enclose plants or inflorescences to stop pollinators visiting the flowers. These physical enclosures could affect estimated pollination benefit by preventing pest access, increasing local air temperature and humidity, or reducing light available to the plant.
Alternative experimental designs that control for the effect of enclosures may involve growing all plants in larger enclosures and introducing managed pollinators or conducting hand pollination treatments Poulsen, 1977). This is unlikely to represent the activity of local pollinator communities and the yield benefits encountered in the field (Woodcock et al., 2019). Different experimental designs therefore represent different compromises between on-farm relevance and the level of experimental control.
These compromises could have far-reaching implications for experimental findings and their application (e.g. see Ainsworth et al., 2008). A greater understanding of how methodological differences translate into pollination dependence estimates could help optimise experimental design and between-study comparisons.
Different crop species vary in their dependence on biotic pollination; some species produce no seeds without external pollination while others show more nuanced responses (Aizen et al., 2009).
Different plant genotypes within crop species also vary in their response to additional pollination (Marini et al., 2015) though this has rarely been considered in economic analyses (but see Fijen et al., 2018 andKlatt et al., 2014). Different plant genotypes are cultivated in different regions, due to a range of factors including consumer preferences and climate requirements. Understanding the range of pollination dependence within a crop (including in breeding material) is necessary to produce a robust estimate of the current and potential future extent of pollination dependence. will increase faba bean yield. Our findings support ecological intensification and specifically the management of pollinators to maximise pollination benefits to faba bean production.

K E Y W O R D S
evidence synthesis, faba bean, pollination dependence, prediction interval, systematic review,

Vicia faba
In this paper, we focus on faba bean Vicia faba L., a partially pollinator-dependent grain legume that has a high protein content and is used as a human food and in livestock feed. Faba bean itself has a valuable role in ecological intensification, supporting the growth of other crops by increasing soil nutrient availability, and supporting pollinating insects by provision of floral resources (Köpke & Nemecek, 2010). The crop produces seed by three mechanisms; autonomous self-fertilisation (autofertility), self-fertilisation mediated by insect visitation and cross-fertilisation following pollen transfer by insects (Drayner, 1959). The relative contribution of each mechanism is known to vary with genotype and the level of cross-fertilisation in previous generations (Frusciante & Monti, 1980). Various authors over the last 60 years have recommended that the crop be supplemented with insect pollinators using domesticated pollinators or pollinator habitat creation (Cunningham & Le Feuvre, 2013;Kendall & Smith, 1975;Riedel & Wort, 1960) but this has not become commonplace.
The variation in pollination benefits to yield, in faba bean and other crops, makes it difficult for practitioners to understand the likely benefits of management for pollinators on a given farm and presents a barrier to ecological intensification and its associated benefits . Meta-analysis is a commonly used approach to synthesise across different sources of evidence and provide an overall effect size estimate (Higgins et al., 2019). However, rather than investigating or accounting for complexity, meta-analyses often simplify data, analysing only a limited subset of available data (or averaging across it) to reduce problems associated with dependence (e.g. repeated measurements across years or different yield metrics; Noble et al., 2017). Recent advancements in meta-analytic approaches now allow us to include more data while accounting for its dependence.
We use such an approach, first using a literature review to identify all experimental estimates of faba bean pollination dependence, followed by application of multi-level meta-analytic models (e.g. Nakagawa & Santos, 2012). We quantify the effect of different factors such as pollinator species on the benefit of pollination via fixed effects. Our approach also allows us to quantify the variation (heterogeneity) associated with other contextual factors via random effects (Nakagawa & Santos, 2012). Furthermore, rather than just simplifying heterogeneity to statistics, we demonstrate that using prediction intervals (IntHout et al., 2016) we can use our understanding of heterogeneity to predict the likely range of pollination dependence across other contexts.

| Identifying publications
To identify relevant publications for our meta-analysis, we searched Web of Science and Scopus using the terms 'faba' and 'pollination' in November 2018. To include grey literature, we conducted a search of Google Scholar using the terms 'faba' and 'pollination' in December 2018; we sorted results by relevance, then checked the first 200 articles (Haddaway et al., 2015). For robustness, we then screened the reference lists of publications that had met our inclusion criteria. We only included publications in our meta-analysis if they compared yield production between plants receiving a pollination treatment and plants that did not receive a pollination treatment. We only included publications that presented yield using at least one of four metrics most likely to correspond to economically relevant yield-yield mass, bean number, pod number or bean number per pod . Several publications reported yield using other metrics that do not translate to economically relevant yield production, these were not included in the meta-analysis but are summarised in Table S1. We included publications that represented insect pollination by hand pollination defined as either tripping (pulling apart keel and wing petals) or hand pollination (tripping + transfer of external pollen). One person (JB) conducted the literature review, we include further details about the review and a PRISMA diagram in the Supporting Information.

| Calculating effect sizes
To quantify the effect of biotic pollination on faba bean yield (the effect size), we used the natural log of the response ratio (lnRR), which is the log proportional change in yield between plants receiving a pollination treatment and undisturbed plants (Hedges et al., 2016).
This converts to the proportion of yield lost without pollination using the simple formula 1 − exp(lnRR).

| Weighting effect sizes
A meta-analysis weights effect sizes by the inverse of their sampling variance. Where standard deviation (or a convertible alternative) was missing from publications, we requested this data from authors for papers dated 2000 onwards. When a measure of variance was not obtainable, we imputed standard deviation based on the fitted relationship between mean, SD and n in the available data (91 data rows imputed; adjusted-R 2 of these models were between 0.29 and 0.46). If a variance measure was not available and we could not impute it, then the effect size was excluded from our meta-analysis (Table 1; Table S1). See Supporting Information for more details and a sensitivity analysis ran without imputed SD; the results were similar.
In medical and eco-evolutionary meta-analyses, n is typically equal to the number of individuals tested in a given experiment. For many effect sizes in our analysis, there was a difference between the number of experimental units (e.g. field patches, or cages containing plants) and the number of individual plants. Where available, we collected information about both sample size types. We conducted our primary analysis using n as the number of experimental units and include findings for n as number of individuals in the Supporting Information; the results were similar. TA B L E 1 Summary of publications reporting effect of pollination treatment on faba bean yield. Values are mean pollination dependence, and in brackets, the minimum and maximum pollination dependence for publications that reported results for different years, genotypes or pollinators. Abbreviations of pollinator types; Apis (honeybee), Bomb. (bumblebee), Open (open pollination), Trip (manual tripping), Hand (hand pollination). Publications marked † did not include appropriate variability information. Publications marked with asterisk were not included in the meta-analysis as variability could not be imputed

| Multi-level meta-analysis models
The publications identified in our literature review were diverse; for example, they were conducted in different countries, they used different pollinator species and different cultivars (plant genotypes), and they measured yield in different ways. The publications often reported more than one effect size; for example, they compared several cultivars, repeated experiments across several years or reported yield using several metrics. The 22 publications that satisfied our inclusion criteria for the meta-analysis (Table 1) included a total of 277 effect sizes.
These differences within-and between-publications allowed us to investigate what causes variation (heterogeneity) in pollination benefit.
Including more than one effect size from the same publication (or from publications that are related in some way, e.g. because they tested the same cultivar) can however lead to correlations (clustering) among these effect sizes, which invalidates model assumptions of independence (Noble et al., 2017). Throughout, we used multilevel meta-analytic models to account for this dependence via random effects and with sampling variance-covariance matrices. These multi-level models follow the same principles as linear mixed effects models (LMMs; Harrison et al., 2018).
To estimate the overall pollination dependence of faba bean, we use a multi-level meta-analytic model with only random effects (a null model, analogous to an LMM with no fixed effects). We identified the optimal random effects structure by comparing AIC of different candidate models. The random effects we tested were an individual effect size identifier (unique per data row, necessary to estimate residual heterogeneity), a publication identifier, year nested within publication (where a publication reported experiments conducted across multiple years), cultivar (plant genotype; a partially crossed-random effect as five cultivars were tested across different publications), the country in which the experimental work took place and the author team (publications sharing an author may have methodological similarities that we did not capture with the other random effects). Apart from the effect size identifier, all the random effects were potential clustering factors. The optimal model is model RE0 in Table 2. We estimated total heterogeneity in pollination benefit and that associated with different levels of clustering (random effects) using I 2 (Higgins & Thompson, 2002).
As well as letting the model estimate dependence of effect sizes due to clustering (random effects, see above), we also explic- Bagging control with 3 levels (controlled, semi-controlled, not controlled)

Fi2
Plant growing conditions with 2 levels (field or plant pot)

Cu2
Type of cultivar with 2 levels (commercial or breeding)

Sc3
Scale of experimental aggregation with 3 levels (within plant, plant, plot)

Ye0
Publication date as a continuous variable TA B L E 2 Summary of all statistical models reported in the manuscript, presented in order of appearance. Likelihood ratio tests were used for all model comparisons. For associated R code, see Table S2 different pollinator species compared to an excluded control group).
We addressed this dependence by dividing sample size for the control group evenly among the shared comparisons before calculating the effect sizes (Higgins et al., 2019). In other cases, multiple effect size estimates were provided per individual (or experimental unit) because multiple yield metrics were measured. We used variancecovariance matrices to account for the resulting correlated sampling (error) variance (Noble et al., 2017). In the main text, we present results from a conservative model that assumed a correlation of 0.5 between effect size sample variances that were measured on the same experimental units. In the Supporting Information, we present results from more conservative models that assume a correlation of 0.9, the results were qualitatively similar.
We added moderators (fixed effects) to the optimal random ef- test. We ran separate models for each moderator and did not test interaction terms because there were many incomplete combinations.
We use marginal R 2 to show how much heterogeneity is explained by moderators (Nakagawa & Schielzeth, 2013).
See the Supporting Information for a detailed description of all tested random effects and moderators, annotated R code, and an extended explanation and justification of all our methods. All model comparisons were made using models fitted with maximum likelihood while results are reported from models fitted with restricted maximum likelihood. All models that we report results from are documented in Table 2, which corresponds to

| Publication bias and sensitivity analysis
Our meta-analysis uses existing literature and it is important to check for biases in that literature (e.g. only publishing significant results). We tested for bias in two ways: (a) visual inspections for a pattern of asymmetry in a funnel plot of the 'meta-analytic residuals' from the meta-regression model including all significant moderators (Nakagawa & Santos, 2012) Table 1 summarises all publications included in our meta-analysis.

| RE SULTS
Several other publications only partially met our inclusion criteria so were not included in the meta-analysis, these are summarised in  Figure 1). If we assume that faba bean crops are optimally pollinated, this indicates that 1.5 million tonnes of faba bean production globally is due to insect pollination (Table 3).
There was very high heterogeneity in pollination dependence both within and between publications (from RE0; I 2 total null model = 0.918, I 2 studyid = 0.054, I 2 yearinstudy = 0.074, I 2 cultivar = 0.484, I 2 infoID (residual) = 0.305). The prediction interval is a highly skewed distribution when transformed to the scale of percentage yield loss without pollination (PI = −69.5% to 73.5%; Figure 1) so we consider it more useful to calculate the probability that a farmer would experience a benefit of biotic pollination to faba bean yield (area under the curve beyond zero; Figure 1), which is 80%.
Plant genotype (cultivar) accounted for more than half of the heterogeneity between effect sizes. If heterogeneity due to cultivars is removed (see Supporting Information), the prediction interval for pollination dependence becomes substantially smaller (−28% to 65%) translating to an 89% chance of seeing a benefit (Figure 1

TA B L E 3
Translating yield dependence into global estimates of pollination value. We calculate this by multiplying the proportion pollination dependence by the global production (4.6 Mt annual mean global production 2009-2018), which comes from FAOSTAT accessed 28/02/2020 The scale of yield measurement did not significantly affect the estimate though there were few effect sizes available to make this comparison, with the majority of effect sizes measured at a wholeplant scale (LRT Sc3 vs. RE0, p = 0.725).
There is weak evidence that pollination benefits vary with pollinator species. While four of the five pollinator groups tested had similar estimates, honeybees were less effective (LRT Po2 vs. RE0, p = 0.110, R 2 marginal = 0.008; Figure 3). We found no difference between estimated pollination depen- Whether the plants tested were commercially available cultivars or breeding lines did not significantly change the dependence estimate (LRT Cu2 vs. RE0, p = 0.618), but there was a greater range of effect sizes in the breeding lines (Figure 4; F test with null hypothesis that variances were equal, p < 0.001, ratio between variances 0.19).
There was little evidence of publication bias. In a model using only non-imputed SD values, and containing moderators of sqrt(vi; where vi = sampling variance for effect size, so sqrt(vi) is standard error), year of publication, and yield metric (3 levels), the slopes of sqrt(vi) and year were not significant (p = 0.791 and p = 0.080 respectively; Supporting Information). Using the dataset that included imputed SD values, there was a significant trend of decreasing pollination dependence over time (LRT Ye0 vs. RE0, p = 0.036; Figure 5).

| D ISCUSS I ON
Our results show that faba bean loses an average of 32.9% yield without biotic pollination. Multi-level meta-analysis allows us to quantify variability (heterogeneity) around the average estimate and attribute this to different sources (Nakagawa & Santos, 2012). The confidence interval, 21%-43%, shows the most likely location of the cross-study average effect and we suggest this information is used in large-scale economic valuations of pollination services (e.g. Table 3).
We found high heterogeneity in pollination dependence. The prediction interval (−70% to 74% in our overall model) incorporates the heterogeneity between effect sizes to illustrate the pollination benefits that a future experiment (or farmer growing a crop of beans) may find (IntHout et al., 2016). Confidence and prediction intervals are distributions, and they become highly skewed on the response ratio (RR) scale (Figure 1). We can, more usefully, calculate that there is an 80% probability that a given farmer will see a benefit of biotic pollination to faba bean yield. Using multi-level meta-analysis and reporting of prediction intervals or their derivatives in this way gives an honest appraisal of biotic crop pollination benefits to farmers.
More than half of the variation in pollination dependence was due to plant genotype. Authors calculating the economic value of biotic pollination commonly account for differences between crop species (Aizen et al., 2009), but there is a clear need to consider differences within crop species (e.g. Klatt et al., 2014) when F I G U R E 4 Orchard plot for commercial versus breeding cultivars. Note that model estimates come from models run on separate subsets of data so that assumptions of homogeneity were not invalidated F I G U R E 5 Time lag bias or a decline in effect size: the year slope (dashed line) along with 95% confidence interval (red dotted lines) and 95% prediction intervals (blue lines) placing economic values on pollination services and/or planning policies and agronomic management. Plant genotype is something that growers can control by choosing to grow a certain cultivar. Several publications in our meta-analysis set out to compare the pollination dependence of different genotypes. Variability in pollination dependence was much greater between (pre-commercial) breeding lines than commercially available cultivars (Figure 4; Supporting Information). If we re-calculate the prediction intervals without the heterogeneity associated with cultivars (e.g. we assume everyone grows a cultivar with average dependence), pollination dependence becomes more predictable (PI −28 to 65%) and the likelihood of a pollination benefit to yield increases from 80% to 89%. Publications have identified faba bean genotypes that are unaffected by pollinator exclusion, and authors have begun to identify heritable traits associated with this autofertility (Torres et al., 1993). Crop breeders might wish to produce faba bean cultivars that do not depend on biotic pollination, as this would remove one possible constraint on production (Marini et al., 2015) but this is difficult to achieve in practice (see Supporting Information). More experimental research is needed to understand how and why pollination dependence varies between plant genotypes and how this interacts with other factors such as maximal yield potential and disease resistance. There has been a reduction in pollination dependence over the 60-year period considered in our analysis, but it is not possible to distinguish between biological changes in dependence and variations in experimental design or publication bias (Koricheva & Kulinskaya, 2019).
The prediction interval overlaps with zero, suggesting that sometimes, biotic pollination will not benefit faba bean yield.
Multiple factors can limit yield simultaneously (see Garibaldi et al., 2018) and experiments, or growers, may find no benefit of pollination because various other contextual factors are limiting yield.
If a plant does not have sufficient resources to mature seeds, then increases in ovule fertilisation may not translate into increases in yield (Garratt et al., 2018). Likewise, if plants are damaged by pests or adverse weather, then this could nullify any benefit of increased ovule fertilisation, or even reduce efficiency if the plant has already invested in a larger number of flowers and seeds (Melathopoulos et al., 2014). Publications repeated across different years specifically discussed variation in weather between years (e.g. Varis & Brax, 1990). Variation in dependence was also reported between soil types (St-Martin & Bommarco, 2016) and following heat stress . We found high heterogeneity between publications, between different years within publications, and high (residual) heterogeneity that we could not explain. This finding means that the average level of pollination benefit will be specific to a grower's context (e.g. their soil type, akin to between publication variation) but other factors will cause additional variation within that context (e.g. differences in weather). Some publications reported negative impacts of pollination, but this was mostly where there was no significant difference between pollination treatments (e.g. Dekhuijzen et al., 1988, but see . These findings suggest that new work to quantify crop pollination dependence (e.g. in other crops or cultivars where meta-analysis is not currently possible due to limited prior research) needs to be repeated across multiple years and/or sites to determine the range of pollination benefit likely to be experienced.
We found weak evidence that the benefit of pollination is greater when the pollinator is not honeybees. This supports ecological intensification, for example, through habitat improvement for wild pollinators, compared to a more conventional input-based approach of introducing honeybee hives on to farmland. Several publications in our analysis directly compared honeybees with other species and found them to be less effective, requiring repeated floral visits to achieve the same podset as some bumblebees (Garratt et al., 2014;Kendall & Smith, 1975, but see Cunningham & Le Feuvre, 2013).
Several authors also reported differences in efficacy between bumblebee species, though we could not test this in our analysis (Kendall & Smith, 1975;Le Guen et al., 1993). Differences in pollinator efficacy likely relate to differences in foraging behaviour, the volume or location of pollen held on the insect, their between-flower movement and the frequency of flower visits made (Marzinzig et al., 2018). Agronomic interventions targeted at specific pollinator species, for instance sowing field margins with plants of a specific floral structure, may be an effective means to attract species of particular benefit to the crop (Garibaldi et al., 2015).
Our results show that there is no yield penalty associated with enclosing plants to exclude pollinators. It is possible that enclosures exerted a combination of positive and negative effects (and this may vary with context, e.g. level of pest pressure), resulting in no difference in overall estimate. We have been unable to investigate more nuanced potential effects of methodologies, for example, non-controlled enclosures were always used with open pollination and commercial genotypes. Regardless, there was no significant difference in pollination dependence between experiments that controlled for bagging effects and those that did not. Given the importance of context and plant genotype, and the clear need to repeat trials across multiple years and/or sites, this finding could help studies to be more efficient (though this should be confirmed in other crop species).
Highly controlled experiments are resource intensive and this limits the number of plant genotypes or contexts that they can cover.
Such experiments may also use pollinators that do not reflect the abundance and activity of insects on-farm. To estimate the benefits of pollination on-farm and to gain useful agronomic insights (e.g. where to establish new pollinator nesting habitat or floral resources), a comparison between open-pollinated plants and enclosed plants is most likely to represent the pollinator community present (e.g. Woodcock et al., 2019). Large-scale, simplified experiments using this design could provide further insights into the role of context on crop pollination benefits. Growers themselves may even be encouraged to conduct these types of experiments  which could increase uptake and suitability of ecological intensification practices .
In summary, we have used a repeatable and rigorous literature review to identify all available comparisons of faba bean yield with and without biotic pollination, and we have used multi-level meta-analysis to quantify and predict faba bean pollination dependence.
Given the importance of context that is increasingly recognised (e.g. Tamburini et al., 2019), quantitative research synthesis is necessary in valuations of crop pollination services. We found high heterogeneity in pollination benefit due to various contextual factors that are both inside and outside of grower control. While the benefits that growers receive from biotic pollination vary across space and time, there is a high likelihood that biotic pollination will increase faba bean yield. This strengthens the case for management of pollinators to maximise those pollination benefits.

ACK N OWLED G EM ENTS
We thank the authors of the original publications used in our research synthesis, particularly those that provided additional variability information. We thank Wolfgang Viechtbauer for support via the metafor documentation, Josh Moatt and Rose O'Dea for making their scripts available online, and Klara Wanelik and the reviewers for helpful comments on the manuscript. S.N. was supported by the ARC Discovery Grant (DP200100367).

AUTH O R S ' CO NTR I B UTI O N S
J.B. conceived the study, conducted the systematic review and metaanalysis and wrote the first draft of the manuscript and Supporting Information documentation. S.N. advised on the analysis, produced the figures and contributed to manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data available via the Centre for Open Science https://doi.org/ 10.17605/ osf.io/b3e9u (Bishop & Nakagawa, 2020).