An examination of ES in real-world systems is needed because controlled experiments that typify most biodiversity-ecosystem function research do not fully represent ecological reality ; specifically, these experiments do not mimic realistic species abundance distributions or species loss scenarios. First, a skewed or ‘hollow curve’ species abundance distribution, meaning that communities are composed of few common but many rare species, is ubiquitous in nature yet biodiversity-ecosystem function experiments have not mirrored this pattern and have instead used sub-stitutive designs that equalize initial abundances among species . Although sub-stitutive designs are arguably preferable for isolating the effects of species richness, they do make experimental communities less realistic and decrease the potential for one dominant species to provide the bulk of ecosystem function. Second, species are lost from communities nonrandomly, with environmentally-sensitive and rare species being at greater risk of extirpation . In contrast, most biodiversity-ecosystem function experiments have assigned species to plots at random to avoid confounding species richness and species composition. It is well-known that, because of this design, comparing high- and low-richness experimental plots reveals the effects of random species loss, which will under- or overestimate the effects of non-random species loss,stackable planters depending on whether species with high or low contributions to function are lost first .
However, biodiversity-ecosystem function researchers have countered that understanding the effects of random species loss is an important starting point, given that future patterns of species loss may be unpredictable . Throughout this paper, we use species loss to indicate local extinctions, because ES are delivered by local populations . Despite an awareness of these systematic differences between biodiversity-ecosystem function experiments and real-world ES, there is no consensus on whether species richness will contribute more or less to function when experimental results are scaled up to larger, more complex systems. Scaling up biodiversity-ecosystem function research is further complicated because analyzing observational data creates challenges not present in experiments. First, few if any real-world ecosystems allow researchers to independently assess services provided by each species in the community. This precludes the use of analytical approaches commonly used in biodiversity ecosystem function research, especially analyses requiring single-species monocultures . Second, communities assemble and disassemble non-randomly with respect to species’ contributions to ES, making it difficult to separate the effect of richness from the effects of species identity. Third, the ubiquity of skewed species abundance distributions in nature makes it hard to separate the effect of abundance from species identity, if the same species are dominant across sites. Because of these issues, no general analytical method exists for biodiversity-ecosystem services studies and, perhaps for this reason, no consensus exists on the importance of different components of biodiversity for ES . Here, we present a novel version of the Price equation that can analyze temporal variance of any ES, so long as the ES can be expressed as a sum of species-level contributions. Our work builds on the original Price equation from evolutionary biology, and its recent adaptations for biodiversity research .
Fox provided the original framework for analyzing temporal variance with the Price equation, and here we extend it so that it can be used with observational data even if species composition is not nested between sites . Our version of the Price equation partitions between-site differences in temporal variance into three additive terms: variance in ES attributable to richness , composition , and context-dependence . Here, we use ‘abundance’ in place of‘context-dependence’ because in our data this term is determined by patterns of abundance fluctuation over time . The richness term establishes an expectation for how species loss and gain would affect ES if species were identical in the ES each provides. The abundance term quantifies the extent to which species present at both sites contribute more to temporal variance of ES at each site. We explore the temporal variance of pollination services using two large, multi-year datasets on pollination provided by bees. Using the new derivation of the Price equation described above, we ask: What is the relative importance of changes in species richness, composition, and abundance to the temporal variance of ecosystem services? Specifically, we compare the relative importance of richness and composition versus abundance. Field surveys. Our study systems consist of the wild bee pollinators of watermelon Matsum. & Nakai and northern high bush blueberry plants, both of which rely on insect pollination for successful fruit production. Over five years , we sampled wild bee communities at 10 commercial watermelon fields in central New Jersey and eastern Pennsylvania, USA. We also, over three years , sampled wild bee communities at 16 commercial blueberry fields in southern New Jersey. In a post-hoc analysis, we confirmed that differences in our results between study systems were not due to length of sampling . Hereafter, we refer to these fields as sites. We ensured that all sites were at least 1 km apart, beyond the typical foraging radius of most bee species in our study .
We did not include the honey bee in our data collection, primarily because in our system the honey bee is a managed species that is kept in hives and moved in and out of crop fields by bee-keepers and farmers. Thus, the temporal and spatial variation in honey bee abundance is driven by hive placement rather than ecological factors. In addition, honeybees are the property of bee-keepers and farmers, so we cannot collect them. Finally, honeybees are present at nearly all sites, so including them would likely increase the relative importance of ‘abundance’,stacking pots making it conservative with respect to our findings to leave the honey bee out of the analysis. Because watermelon is an annual species, farmers do not necessarily plant it in the same locations each year. To exclude potential effects of spatial variation on wild bee communities, we included watermelon sites in our analyses only if the maximum among-year distance between transects was ≤ 435 meters; this is within the typical foraging radius of all but the smallest bee species in our study . Pollination services. To measure bee richness and the pollination services delivered at each site on each date, we collected two forms of data: the number of individual wild bees visiting flowers, and the number of pollen grains deposited per flower visit. We then multiplied each species’ abundance by the mean pollen deposition of its morphological group to obtain that species’ contribution to pollen deposition. To measure bee abundance, we established a 50 m transect at each watermelon or blueberry site. We collected bees visiting flowers by net throughout the transect and then processed voucher specimens for species-level identification by taxonomists . At each site, data were collected on three days during each plant species’ peak bloom period, with three temporally-stratified 20-minute-long collections during the day; all data collection days were sunny, partly cloudy, or bright overcast with limited wind. Total collection effort was 135 days for watermelon and 144 days for blueberry. We measured the pollination efficiency per flower visit for different pollinator groups in field experiments. We offered a virgin flower to an individual bee foraging in the field, allowed the bee to visit the flower one time, and recorded the pollinator group of the bee . These pollen deposition experiments were conducted in three years for watermelon, and in two years for blueberry. Back in the laboratory, we use a compound microscope to count the number of con-specific pollen grains deposited during the single flower visit. To prepare slides, watermelon stigmas were softened in 10% KOH, and stained with 1% fuchsin. Blueberry stigmas were softened in 1 M NaOH, and then stained for 48h in 0.01% analine blue buffered in 1 M K3PO4. In the pollen deposition experiment, individuals were placed into morphologically similar groups which could be differentiated in the field . Due to the difficulty of collecting single-visit pollen deposition data, we did not collect these at all sites and in all years. However, we do know that that different morphological groups differed significantly in pollen deposition rates, while species within groups did not . Further details of all data collection methods, and site details, are available for both watermelon and blueberry . Because there is substantial variability within each morphological group’s distribution of pollen deposition rates, we conducted a sensitivity analysis to test whether choosing a single pollen deposition value and discarding the remaining variability affected our results. Instead of multiplying each bee species’ abundance by its morphological group’s mean pollination efficiency, for each individual bee, we randomly drew one pollination efficiency value from the correct morphological group’s distribution.
We repeated this sensitivity analysis 1000 times, and found that the same results were obtained using either the ‘mean’ or ‘sensitivity’ versions of the analysis . Temporal variance of pollination services. Our data exists along three axes: sites, species, and years. We organized this data in a series of matrices, with one matrix for each site. Specifically, let Si be an m x n matrix describing the contribution of bee species n to pollen deposition in year m at site i. The contribution of bee species n equals its abundance times its group’s mean pollination efficiency . Then, let Vi be the variance of Si, which is an n x n variance-covariance matrix, with the diagonal elements describing the variances of each species, and the off-diagonal elements representing the covariances between any pair of species. The sum of Vi is then the total temporal variance of ES at site i. The Price equation decomposes the differences in total temporal variance of ES between any two sites, e.g. between V1 and V2. Total temporal variance must be non-negative, but an individual species can have negative contributions to total temporal variance if its summed covariance is negative and greater in magnitude than its variance. Thus, there is no restriction on the shape of the relationship between species richness and total temporal variance, except that total temporal variance must be non-negative. Although most biodiversity-stability studies have measured stability using the inverse of the coefficient of variation , the CV is not always a good indicator of stability because it conflates how mean ES and temporal variance of ES respond to the same environmental changes . In addition, the CV of ES cannot be expressed as a sum of independent species-level contributions and thus a broadly applicable Price equation partition is not possible, although versions can be developed for particular cases . Price equation partition. Our partition combines the advances of two previous versions of the Price equation. First, Fox showed the general framework for using the Price equation to analyze temporal variance. Second, Fox and Kerr , presenting a Price equation partition of mean ES , showed how to relax a previously-existing ‘nestedness’ requirement. In ecology, the nestedness requirement means that, when comparing two sites, only one site could contain a species not found at the other site . That step is crucial for real-world applications because species composition is rarely nested. Here, we apply the logic of Fox and Kerr —specifically, how to remove the nestedness requirement—to extend Fox’s method for analyzing temporal variance with the Price equation. The result is the first version of the Price equation suitable for analyzing temporal variance of ES in real-world systems. In the watermelon system, we collected 3044 individual wild bees belonging to 59 species , 14 genera, and 5 families. Sites ranged in total temporal variance from 2.3 million to 51.1 million grains of pollen , in richness from 20 to 32 species , and in abundance from 140 to 435 individuals ; richness and abundance values are summed across years. In the blueberry system, we collected 1067 individual wild bees belonging to 36 species , 9 genera, and 4 families. Sites ranged in total temporal variance from 2.6 to 875,154 grains of pollen , in richness from 3 to 19 species and in abundance from 4 to 143 individuals . The range of richness values is similar to the range found in experiments that study various ecosystem services , although only two of our sites had species richness values similar to the lower end of richness manipulations . The rank-abundance distributions for each system are shown in Figure 1. The five Price equation terms are measures of directional effect size, where positive or negative values tell how strongly each term decreases or increases between-site differences in the total temporal variance of ES. For both systems, all richness-loss terms were negative and all richness-gain terms positive, with richness-loss greater in magnitude .