Engineering a plant community to deliver multiple ecosystem services

. The sustainable delivery of multiple ecosystem services requires the management of functionally diverse biological communities. In an agricultural context, an emphasis on food production has often led to a loss of biodiversity to the detriment of other ecosystem services such as the maintenance of soil health and pest regulation. In scenarios where multiple species can be grown together, it may be possible to better balance environmental and agronomic services through the targeted selection of companion species. We used the case study of legume-based cover crops to engineer a plant community that delivered the optimal balance of six ecosystem services: early productivity, regrowth following mowing, weed suppression, support of invertebrates, soil fertility building (measured as yield of following crop), and conservation of nutrients in the soil. An experimental species pool of 12 cultivated legume species was screened for a range of functional traits and ecosystem services at ﬁve sites across a geographical gradient in the United Kingdom. All possible species combinations were then analyzed, using a process-based model of plant competition, to identify the community that delivered the best balance of services at each site. In our system, low to intermediate levels of species richness (one to four species) that exploited functional contrasts in growth habit and phenology were identiﬁed as being optimal. The optimal solution was determined largely by the number of species and functional diversity represented by the starting species pool, emphasizing the importance of the initial selection of species for the screening experiments. The approach of using relationships between functional traits and ecosystem services to design multifunctional biological communities has the potential to inform the design of agricultural systems that better balance agronomic and environmental services and meet the current objective of European agricultural policy to maintain viable food production in the context of the sustainable management of natural resources.


INTRODUCTION
Sustainable socioecological systems rely on the integration of multiple ecosystem services at the scale relevant to the underlying ecological processes (Kremen 2005, Bennett et al. 2009, Carpenter et al. 2009, Diaz et al. 2011. Focusing exclusively on a single service risks instability and loss of function. Nowhere is this more apparent than in agriculture, where the intensification of food, energy, and fiber production (provisioning services) has been pursued at the expense of regulating and supporting services such as pollination, bio-control, healthy soil, and clean water (Power 2010, Raudsepp-Hearne et al. 2010. Agricultural policy in Europe is currently attempting to redress this balance by structuring agricultural subsidies in a way that explicitly links food production to the maintenance of ecosystem services (Mouysset 2014). Agriculture in the European Union is supported by the Common Agricultural Policy (CAP) via two ''pillars'': pillar 1 for supporting food production, and pillar 2 for rural development (including voluntary agri-environment measures). In the latest reform of the CAP in 2014, support for production under pillar 1 includes a compulsory ''greening'' element for the first time; a proportion of the cultivated land is required to be managed for ecosystem services, and rotations need to include a minimum diversity of crops (European Commission 2013). However, in the context of a growing world population, managing land for ecosystem services needs to be done in such a way that Manuscript  food production is not compromised, and methodologies need to be developed that allow the delivery of sometimes conflicting ecosystem services to be quantified and reconciled (Nelson et al. 2009). One approach to quantifying ecosystem service delivery, for which a number of case studies now exist, is to use metrics of functional traits to predict the ecosystem function of different communities (Diaz et al. 2007, de Bello et al. 2010, Lavorel et al. 2011. The usefulness of the approach for modeling multiple ecosystem services delivered by plant communities, including productivity and soil nutrient cycling, is being increasingly demonstrated (Minden 2011, Pakeman 2011, Laliberte and Tylianakis 2012, Lienin and Kleyer 2012, and the conceptual framework has the potential to incorporate services delivered by higher trophic groups . In this study, instead of using these models to quantify the functionality of existing semi-natural habitats, we use them to inform the design of a cultivated plant community. The relationship between species richness and the delivery of a single ecosystem service is case specific and not always positive (Hooper et al. 2005, Balvanera et al. 2006. However, when multiple services are assessed in parallel, increasing species richness may be seen as desirable insofar as different species perform complementary functions (Hector and Bagchi 2007, Gamfeldt et al. 2008, Zavaleta et al. 2010. Where overlap between species in terms of their contribution to different services is small, the multifunctionality of the system has been predicted to continue to increase as additional species are added to the community (Hector and Bagchi 2007). However, where there are trade-offs between services or where species contribute negatively to a service, the relationship between species richness and multifunctionality may quickly saturate or become negative (Raudsepp-Hearne et al. 2010, Zavaleta et al. 2010, Gamfeldt et al. 2013). If it is also assumed that the level of any given service is also largely determined by the attributes of the dominant species in the community (Grime 1998), multifunctionality will not only be determined by the combination of services and functional diversity of the species pool, but also by the dominance hierarchy and competitive dynamics of the community. We took a process-based approach to modeling these competitive interactions, combined with functions quantifying the relationships between functional traits and ecosystem services, to engineer a community of cultivated legume species that delivered a balance of six ecosystem services: early productivity, regrowth following mowing, weed suppression, support of invertebrates, soil fertility building (measured as yield of following crop), and conservation of nutrients in the soil.
Legume-based cover crops are currently managed to deliver high levels of biomass of forage with high digestibility and to build soil fertility. However, the current reliance on a few highly productive legume species results in residues that are rapidly mineralized after ploughing and an asynchrony between nutrient supply and demand from the following grain crop, with consequent losses of nitrogen to the atmosphere and water courses (Crews and Peoples 2005). There is, therefore, potential to design species mixtures for cover crops that better reconcile agronomic and environmental functions. A range of legume species were grown in monocultures at multiple sites across a geographical gradient in the United Kingdom in order to screen for functional traits and ecosystem service delivery. These data were used to identify the species mixture with the appropriate level of complexity that reconciles productivity with a more recalcitrant residue composition and a number of other ecosystem services: soil fertility building (assessed as following grain crop yield), weed suppression, and support of invertebrates. To validate our approach, a complex mix of ten legumes and four grass species, referred to as the ''all species mix'' (ASM), was also grown at all sites and assessed for the delivery of the ecosystem services. The specific system analyzed here is presented as proof of concept, but our approach is relevant to any cultivated species mixture, including pastures (Finn et al. 2013), intercropping (Damour et al. 2014), and agri-environment habitats (Balzan et al. 2014).

MATERIALS AND METHODS
Our methodology followed a number of logical steps.
(1) We chose candidate species based on expert knowledge and trait databases. (2) We grew monocultures of the chosen species in the field at multiple sites to quantify ecosystem services and functional traits. (3) We simulated plant growth and competition for all potential mixtures of species. (4) We identified the optimum number and combination of species using an index of multifunctionality. (5) Finally, we validated predictions from the competition model using data from plots sown with a complex species mix at all sites.

Choice of candidate species
In seeking to optimize delivery of multiple ecosystem services, it was desirable to include in the field experiments species that represented the range of functional space occupied by the available legume flora. An initial list of 22 candidate legume species was compiled, and data on biological and agronomic variables was obtained from the literature and expert knowledge (Table A1). A principal components analysis was done to identify functionally dissimilar species. In addition, the tolerance of each species to grazing, autumn sowing, and frost damage was used as additional agronomic filters; any species that was intolerant of any two factors was excluded. Twelve species from across the ordination space were then screened in the field for the delivery of multiple ecosystem services. Legume-based cover crops are grown for a number of agronomic services, including biomass production for forage or green manure, soil fertility building, and weed suppression. In addition to the measurement of these services, two environmental services were also measured in the field: support of invertebrates, and reduction of nitrogen leaching through a more recalcitrant residue composition. Published literature on the functional traits related to these services was used to inform the selection of traits measured on the experiments.

Field experiments to quantify ecosystem services and traits
The legumes were sown in April 2009 in monocultures in plots with a minimum size of 5 3 1.2 m in a fully randomized block design with three replicates. The experiment was repeated at five field sites with a wide geographical coverage across the United Kingdom: Rothamsted Research, Hertfordshire (51848 0 38 00 N, 0822 0 02 00 W); Duchy College, Cornwall (50813 0 38 00 N, 5818 0 23 00 W); Wakelyns Agroforestry, Suffolk (52821 0 37 00 N, 1821 0 09 00 W); the Scottish Agricultural College (SAC), Aberdeen (57811 0 06 00 N, 2812 0 45 00 W); and Aberystwyth University, Wales (52825 0 48 00 N, 4801 0 22 00 W). Time of emergence and final plant density were assessed on two 0.25-m 2 fixed quadrats on each plot. The legumes were mown in early summer and autumn and either incorporated in the autumn of 2010 or spring of 2011. Lathyrus pratensis established poorly at a number of sites, and Vicia sativa was killed by mowing; they were, therefore, excluded from further analysis.
Six ecosystem services were assessed on the experiments. Early productivity, regrowth following mowing, and weed suppression were measured at all sites. Productivity was measured both as early biomass and regrowth, as the aim was to optimize provision of forage at both the first and later cuts. Early productivity was assessed before the first mowing in July/August by fitting an exponential function to a time series of individual plant dry mass sampled at five intervals before mowing and converted to grams per square meter on a standard date, 15 July, using plant density data. After mowing, the aboveground biomass was sampled at weekly intervals from five separate 0.25-m 2 quadrats and expressed as specific aboveground net primary productivity (SANPP, gÁg À1 Ád À1 ; Vile et al. 2006). The final regrowth sample was also used to assess the relative ability of the legume species to suppress weeds by separating out the weed biomass and measuring weed dry mass.
Three further services were measured at one site, Rothamsted (Hertfordshire, UK): support of invertebrates, which is important as a reservoir of natural enemies of crop pests and a food resource for farmland birds (Storkey et al. 2013); soil fertility building; and the conservation of nutrients in the soil. Invertebrates were measured using a vortis suction sampler (Arnold 1994). Samples consisted of five 10-second ''sucks'' taken from the legume monoculture plots, and total invertebrates were counted in each sample. Next, soil fertility building was assessed indirectly by measuring the yield of a winter wheat crop, grown in the absence of additional fertilizers, following the incorporation of the legume residues in September 2010. Yield was measured on a 1m 2 quadrat sampled by hand from each plot. The final service of interest was the conservation of nutrients within the soil, mitigating against diffuse nitrogen pollution. It was not possible to measure the breakdown characteristics of the legume residues directly in the context of reducing nutrient losses from the system. However, there is an established relationship between the (lignin þ polyphenol) : N ratio of legume residues and the rate of N mineralization (Fox et al. 1990). Sufficient biomass remained for all species immediately before incorporation to analyze the lignin, polyphenol, and nitrogen contents of the residues at Rothamsted.

Simulation model of plant growth and competition
A simulation model of growth and competition of multiple plant species has previously been developed to predict competition in wheat from annual weeds (Kropff andSpitters 1992, Storkey andCussans 2007). The model, written in Cþþ, was parameterized for the 10 legume species and adapted to simulate regrowth of the canopy after mowing. The model simulates growth on a daily time step and was divided into three phases. Before the onset of competition for resources, plants are assumed to be sink-limited and growing exponentially according to a relationship between relative growth rate (RGR) and thermal time where W is aboveground dry mass, W 0 is initial mass, T is daily mean temperature, and T b is base temperature. The exponential growth phase was parameterized for total dry mass, aboveground dry mass, and green leaf area by sequentially sampling seedlings grown in pots in April/May 2010 and 2011, according to previously published protocols (Storkey 2004). The model assumes a total green area index of 0.75, representing the onset of competition for resources and the start of the second phase in the model (Kropff and Spitters 1992). After this point, growth was modeled from the radiation intercepted by each species, assimilation rates, and conversion efficiencies. The model integrates the radiation intercepted by the competing species in each of five horizontal layers in the canopy from the species-specific light extinction coefficients and vertical leaf area distribution, where I a,h,i is the light absorbed by species i (JÁm À2 Ás À1 ) at height h (m), q is the reflection coefficient of the canopy, k i is the extinction coefficient for species i, I 0 is incident radiation (JÁm À2 Ás À1 ), k j is the extinction coefficient of species j, and L h, j is the leaf area index of species j at height h ( j ¼ 1, . . . , n species in the mixed canopy [Kropff 1993]). Radiation intercepted at each layer in the canopy was used to calculate the instantaneous assimilation rate (kg CO 2 Áha À1 Ás À1 ) based on initial light use efficiency and maximum assimilation rate, A max , using a generic function that relates A max to leaf temperature, leaf N concentration, and specific leaf area (SLA, g/m 2 ; Storkey 2005). Assimilate is converted to biomass of different plant organs in the model, using conversion efficiencies based on C:N ratios (Penning de Vries et al. 1974) and partitioning functions plotted against photothermal time. The second phase of the model required parameters for height growth, partitioning, vertical leaf area distribution, and C:N ratios of different plant organs. All these parameters were measured on the monoculture plots at the Rothamsted site, using established screening protocols (Storkey and Cussans 2007).
The parameters for the early exponential growth phase before the onset of resource competition will not be affected by the species identity of the in silico plant mixtures used in the optimization exercise. However, it is likely that partitioning parameters later in the season may differ according to the identity of the neighboring plants despite the fact that the plant densities of the modeled mixtures are equivalent to the experimental monoculture plots used in the parameterization. Given the complexity of modeling this phenotypic plasticity, however, and for the purposes of the optimization exercise, we assumed the phenotypic response to interspecific competition was equivalent to the observed response to intra-specific competition.
The third phase of the model simulated regrowth after mowing in the summer or autumn. Loss of biomass as a result of mowing was calculated as a function of plant height, vertical distribution of foliage, and mowing height. Detailed physiological measurements were not taken in the field to parameterize subsequent regrowth; instead, sequential biomass samples were taken postmowing in the monoculture plots and plotted against accumulated radiation intercepted, calculated in order to quantify radiation use efficiency for each species. Biomass production post-mowing in the mixtures was modeled using the function for light interception of each species and the species-specific value for radiation use efficiency. A separate function describing height postmowing was used to model competition for light.
The parameterization of the eco-physiological model of competition generated values for a range of plant functional traits, namely, seed mass, maximum height (measured at all sites), specific leaf area (SLA; g/m 2 ), leaf : stem ratio (L/S), C:N ratio of mature leaves, and stems and leaf N content. These trait data, along with the residue composition data, were analyzed to quantify relationships between legume traits and ecosystem services for 10 species, for which data were available from all sites. For each service, all subsets linear regression was used to identify the combination of independent traits that explained the maximum variability using only explanatory variables with P , 0.05. Where a service was measured at multiple sites (early productivity, regrowth, and weed suppression), the analysis was done on the mean value. As opposed to a step-wise approach, all subsets regression analyzes all possible combinations of explanatory variables, using the adjusted R 2 and Mallows C p as criteria for comparing models and inclusion of covariates. However, the assumption of linearity of the relationships of traits to services is a potential limitation of this approach.

Identification of optimum number and combination of species
The competition model predicted relative biomass of the component species in all possible species mixtures. This output was used to calculate the community weighted mean (CWM) of each of the ecosystem services, measured at all sites using site specific data on service delivery in the monocultures. For the remaining services, measured at Rothamsted only, the models of the relationship with functional traits derived above were used to predict services at the other sites, using CWMs calculated using the competition model output. The exception was N mineralization rate, which was not measured directly in the experiments; the CWM of the (lignin þ polyphenol) : N ratio for the residues of the different mixtures was, therefore, used in the analysis as a proxy for this service. To standardize the data, the relative performance of the mixtures was expressed as a proportion of the best performing mix for each service (with a maximum of one) at each site.
It was assumed that all services were equally important and that any index that combined their contributions had to be limited by the lowest level service. The combination of multiple services, therefore, was considered using a limiting factor approach, which has previously been modeled using a sum of reciprocals function (Aikman and Scaife 1993: Eq. 3). This is a way of combining multiple limitations such that the result cannot exceed the smallest component while reflecting the effect of all constraints; if a particular mix was the best performer for all six services, this would result in a maximum value for I of 0.167 where S i is the relative performance of a mixture for service i, compared to the highest value for that service.
It is possible to identify the optimum number and combination of species from the community with the maximum value for I. However, when interpreting the results, the sampling effect must be taken into account; there are many more combinations of species at intermediate levels of species richness than at high or low levels, increasing the probability of deriving higher maximum values of I at intermediate species richness. To separate the underlying biological processes from this statistical artefact, two further analyses were done. First, a null model was run in which the 2 10 À 1 ¼ 1023 possible combinations were shared equally between the 10 levels of species richness, using dummy species with trait values sampled randomly from frequency distributions fitted to the observed data from the experimental species pool. For each simulated community, average trait values were used to calculate the delivery of each service, and a constant biomass was used in the calculation of following crop yield and weed suppression. Second, the null model was rerun using the same, unequal numbers of communities at each level of diversity, as in the analysis of the experimental species pool (10,45,120,210,252,210,120,45,10, and 1 for 1-10 species in the community, respectively).

Validation using data from all species mix
To validate our approach of using the simulation model to predict trait metrics and services delivered by mixtures, an additional plot containing 10 of the legume species, defined as the all-species mix (ASM), was sown in the field experiments at all sites (see Plate 1). The ASM also contained four grass species (Lolium multiflorum, L. perenne, Phleum pratense, and Festuca pratensis) to reflect the common practice of including a proportion of grass in legume cover crops to enhance nitrogen fixation. Relative abundance of each species in the ASM at each site was calculated from the proportion in the seed mix adjusted by site-specific emergence counts from the monocultures that were also used to parameterize relative time of emergence. The relative biomass predicted by the competition model at different points in the season was then used to predict the CWM Notes: Early productivity, regrowth post-mowing, and weed biomass were measured at all sites, numbers of invertebrates and yield of following crop only at Rothamsted site. Abbreviations are L (leaf ), S (stem), SLA (specific leaf area (g/m 2 ). Data from Ibers site is not included in the regrowth analysis because of poor regrowth on all plots. PLATE 1. Example of the all-species mix with ten legume and four grass species sown at all field sites to validate the simulation model. Photo credit: J. Baddeley. of functional traits and ecosystem service delivery from the regression models, in the same way as the optimization, and compared to observed data.

RESULTS
Highly significant relationships were found between ecosystem services and functional traits in the legume system (Table 1). There were also trade-offs between functional traits and, consequently, between the ecosystem services that each legume species delivered (Fig. 1). The index of multifunctionality identified two species, Medicago lupulina (L.) and Trifolium pratense (L.), that were good all-rounders, one of which always featured in the optimum mix at all sites (Table 2). Other species were positioned toward the extremes of the ordination space and performed well on some services but poorly on others, lowering their multifunctionality index. Across the five sites, a relatively  simple mixture of one to four species resulted in the optimum balance between the services ( Table 2). The fact that the highest value of I was observed at intermediate levels of species richness was shown to be partly a result of the sampling effect (Fig. 2b, c). In the null models, with a large species pool distributed evenly over the trait space, there was a theoretical optimal combination of traits for achieving the best balance of services. When the sampling intensity was constant across all levels of species richness, the model could be optimized using a single species with the ideal combination of traits (Fig. 2b). However, it is unlikely that, given trade-offs between traits as shown in Fig. 1, the ideal species will be present in nature. Using the unequal sampling effort that reflected the small species pool used in the field experiments, the likelihood of approaching the optimal combination of traits increased with the number of sampling events (Fig. 2c).
In the field experiments, the diversity of the ASM plots tended to decline with time at all sites, with four species becoming dominant: Trifolium repens, T. pratense, Medicago lupulina, and M. sativa. This was reflected in the output of the eco-physiological simulation model, which predicted that the proportion of the remaining species would decline over time (Appendix: Figs. A1, A2). When the model output of relative biomass of the species in the ASM was used to generate values for the CWM of functional traits and combined with the regression models (Table 1), the delivery of the multiple services across the sites by the species mixture was successfully predicted for most services (Fig. 3). This supports our general approach of combining trait/ service relationships with the competition model, despite the fact that the simulation model could not be parameterized for multispecies mixtures. However, the model underestimated regrowth and weed suppression in the ASM; one possible explanation may be error associated with the simulated growth of the grasses. The competition model was parameterized in monocultures without additional nitrogen, but the observed growth rate of the grasses was higher in the ASM because of facilitation from the legumes. Addressing this problem in future versions of the model will also allow the optimum grass/legume ratio in a mixture to be quantified. The optimum legume species mix (derived from the average delivery of each service across the five sites) was predicted to outperform the ASM for four of the five services measured directly in the field.

DISCUSSION
Delivering multiple ecosystem services from the same plant community will depend on exploiting the functional contrasts between species that both enable them to coexist within the same ecological niche and to deliver complementary ecosystem services. The optimum solution in terms of the number and combination of species will be unique to the functional composition of the   FIG. 2. (a) Index of multifunctionality calculated as the sum of the reciprocals of each service (averaged over the five sites) for all 1023 possible combinations of the 10 legume species (with a maximum of 0.167). Services included in the analysis were early productivity, regrowth after mowing, weed suppression, soil fertility building (measured as yield of following crop), support of invertebrates, and soil nutrient retention; (lignin þ polyphenol) : N ratio of residues was used as a proxy for this last service. Dashed line is the mean of the index. (b) Output of null model with equal number of communities sampled at each diversity level and using dummy species with traits sampled from frequency distributions fitted to observed data from experimental species pool. (c) Null model using unequal numbers of communities sampled at each diversity level, as in (a). Points have been jittered along x-axis. available species pool, the competitive dynamics of the community, and the specific combination of services. In cultivated species mixtures, there is a further constraint of availability of candidate species for inclusion in a seed mix, and pragmatic solutions are required that balance ecological principles with management constraints. Where the available species pool is limited, the optimal level of diversity will be determined largely by the number of candidate species included and the functional space they occupy, emphasizing the importance of the initial step of selecting the candidate species for calculating I. In our system, screening a wider pool of legumes could potentially have identified a single species with an optimal combination of traits or a better performing mix. The development of global trait databases (Kattge et al. 2011) and a growing literature on relationships between traits and multiple ecosystem services (de Bello et al. 2010) is making this initial screen of candidate species a realistic possibility for a range of systems.
For the experimental pool of 10 cultivated legume species used in our study, the optimal solution exploited the temporal and spatial contrasts in the growth pattern of the candidate species in order to optimize the delivery of different services. The best mixes always included a species from the center of the ordination space, complemented, at some sites, with additional species that exhibit contrasting traits, exploiting differences in growth habit. For example, the vigorous growth post-mowing of Medicago lupulina FIG. 3. Ecosystem services were measured on the all-species mix (ASM) plots to validate our approach of predicting ecosystem function in mixtures. Observed values for five services measured at (a, b, e) all five sites or (c, d) only at Rothamsted were plotted against predicted values from relationships with functional traits (Table 1). (a) Early productivity (g/m 2 ), (b) regrowth (gÁg À1 Ád À1 ), (c) yield of following crop (Mg/ha with 85% dry matter), (d) Numbers of invertebrates (no./m 2 ), and (e) weed biomass (g/m 2 ). Open circles are mean values from monoculture plots; solid circles are mean values from ASM plots; predicted values were calculated using proportional biomass output from the simulation model (using site-specific management and weather data as inputs) to calculate CWM (community weighted mean) of the traits used in the regression models in Table 1; dashed lines show predicted level of each service for optimal mixture based on ecosystem service data averaged across sites.
complemented the early productivity of Medicago sativa or Trifolium incarnatum. There is scope within the model for adjusting relative densities of the component species in the mix to further refine performance, and the use of the sum of reciprocals function for combining services allows individual services to be weighted. If, for example, the mix was only going to be used as a short-term ley, more weighting could be put on early productivity. An important service that was not included in our analysis was the support of pollinator communities, for which there were insufficient measurements made on relevant flower traits. Although the optimum mix represented a range of flowering times, it is also likely that a diversity of flower architecture will also be important for supporting a range of different pollinator groups. In this regard, functional divergence may be a more appropriate metric (de Bello et al. 2010).
The specific case study developed here is relevant to farming systems that incorporate legume-based cover crops into the rotation. These are currently dominated by low-input systems, but increasingly conventional operations are considering cover crops in response to increasing weed pressure and cost of inputs. However, the approach we have developed has general relevance to any cultivated multispecies community where there is a need to reconcile production with supporting and regulating ecosystem services, currently an important driver of European Agricultural Policy (Mouysset 2014). In mixed farming systems, short-term pastures are currently managed almost exclusively for productivity, with the emphasis on fast-growing grasses such as Lolium sp. (rye grass). A more functionally rich grassland seed mix, engineered using the concepts we develop in our study, could potentially mitigate some of the environmental problems associated with simple rye grass swards, including increasing carbon storage and pollen and nectar resources and reducing greenhouse gas emissions (Pilgrim et al. 2010). Similarly, complex seed mixes are currently sold to be sown as part of agrienvironment schemes designed to deliver environmental benefits from areas of uncropped land on farms (Balzan et al. 2014). These seed mixes tend to be targeted at single ecosystem services and are marketed as such; examples include wild bird seed mixes or pollen and nectar mixes. There is increasing pressure on farmland to secure food production, and there is, therefore, a strong driver to minimize the amount of land that needs to be taken out of production in order to maintain important ecosystem services delivered by biodiversity. This highlights the potential to apply our framework in order to engineer seed mixes that optimize several of these services from the same plant community (Holland et al. 2014).