Diversifying bioenergy crops increases yield and yield stability by reducing weed abundance.

Relationships between species diversity, productivity, temporal stability of productivity, and plant invasion have been well documented in grasslands, and these relationships could translate to improved agricultural sustainability. However, few studies have explored these relationships in agricultural contexts where fertility and weeds are managed. Using 7 years of biomass yield and species composition data from 12 species mixture treatments varying in native species diversity, we found that species richness increased yield and interannual yield stability by reducing weed abundance. Stability was driven by yield as opposed to temporal variability of yield. Nitrogen fertilization increased yield but at the expense of yield stability. We show how relationships between diversity, species asynchrony, invasion, productivity, and stability observed in natural grasslands can extend into managed agricultural systems. Increasing bioenergy crop diversity can improve farmer economics via increased yield, reduced yield variability, and reduced inputs for weed control, thus promoting perennial vegetation on agricultural lands.

INTRODUCTION

Perennial bioenergy crops provide additional environmental benefits compared with annual bioenergy crops, including increased soil carbon sequestration (, ), reduced nitrous oxide emissions (), less nitrate leaching to groundwater (), and mediated tile drainage flows () through year-long plant cover. However, lower biomass yields of perennial versus annual bioenergy crops can create economic barriers to grower adoption (). Economic viability of perennial bioenergy crops may be achievable through diversification strategies that increase yields, stabilize yields across years, and lower input costs compared to annual bioenergy crop monocultures. Identifying agronomic practices that increase yield and yield stability of perennial bioenergy crops is critical to their economic sustainability, potential for expanded acreage, and thus realization of their environmental benefits.

Many ecological experiments have shown that increasing plant species diversity results in increased productivity (, ). Results from these experiments have been used to promote increasing species diversity of perennial bioenergy cropping systems. However, it is uncertain if and how diversity-productivity relationships observed in ecological experiments translate to agronomic situations because ecological experiments (i) often maintain plant diversity by removing volunteer species and invading weeds, (ii) are not harvested annually using schedules for bioenergy production, and (iii) are not managed to maintain soil fertility in response to nutrient removal during annual biomass harvest (). Relatively few studies have tested the diversity-productivity relationship while managing species composition and fertility for agricultural production [e.g., (–)]; however, these studies indicate that such relationships can increase yields at the field level and food security at a systems level ().

The temporal stability of productivity, which can be defined as the mean productivity through time divided by the SD of productivity through time (), is often increased in communities with higher plant species diversity. This metric is analogous to yield stability in a perennial cropping system. If positive diversity-stability relationships extend to bioenergy cropping scenarios, increasing diversity of bioenergy crops could function as insurance for protecting growers from yield fluctuations during atypical growing seasons (). However, the mechanisms underlying the diversity-stability relationship are complex and not totally clear (), and it is uncertain whether the positive relationship persists under the agronomic conditions of bioenergy cropping systems (, ).

Species asynchrony—the change in relative abundance of species within a community in response to changing environmental conditions through time—has been identified as an important contributor to the diversity-stability relationship (, ). High species asynchrony can occur in a diverse community because (i) individual species varying in traits will inherently respond differently to environmental changes through time and/or (ii) competition among individual species in a community can alter relative abundances, with some species increasing and others consequently decreasing (). Although many ecological studies have consistently identified asynchrony as a positive driver of stability (–), the potential role that asynchrony can play in driving the diversity-stability relationship remains unclear.

Priorities for weed management differ for perennial bioenergy cropping systems compared with annual cropping systems and perennial forage systems, and the influence of plant community variables like species diversity on weed abundance is poorly understood. Annual cropping systems can use weed management strategies including crop rotations to interrupt patterns of common weed growth and reproduction, tillage to physically destroy weeds before or after planting, or herbicides that target weeds with growth forms that differ from the main crop; all of which may not be practical or as useful for mixed-species perennial bioenergy crops. Unlike perennial forage systems, perennial bioenergy systems may not receive multiple defoliation events within a single growing season, which can affect weed population dynamics and reduce the weed seedbank (). Increasing diversity of species mixtures that promote asynchrony can limit available niche space for weeds (, ). We hypothesized that more diverse, species-rich bioenergy cropping systems would result in higher asynchrony of species through time, which would be associated with lower weed abundance and greater yield and yield stability.

We conducted a dedicated bioenergy cropping systems experiment at nine locations varying in weather and soil type (Fig. 1) for 6 or 7 years that included 12 mixtures of native perennial species varying in planted species richness. This study mimicked agricultural conditions in that the species were selected for yield and complementarity, grown with and without nitrogen fertilizer, harvested annually after plant senescence as recommended for energy conversion via gasification or fermentation, and managed with judicious weed control to prevent perennial weed infestations only during the first 2 years and allowed successional trends in planted species to occur with weeds thereafter. By documenting the relative abundance of all species in each experimental unit annually, we used observed species richness, weed abundance, and species asynchrony to explain variation in yield, temporal variability of yield, and stability.

Fig. 1.

Variation in monthly temperature, precipitation, and soil type.

Violin plots showing variability in mean daily temperature by month (A) and mean total monthly precipitation by month (B) across nine experimental locations and 7 years. Soil classification of each experimental site by sand, silt, and clay fraction (C). Points on the soil texture classification chart are proportional to the crop productivity index, a measure of soil fertility ranging from 0 to 100 (see Materials and Methods for details). Additional geographic and soil information can be found in table S1.

RESULTS AND DISCUSSION

Planted species richness in mixed-species perennial bioenergy cropping systems directly increased yield, reduced weed abundance, increased asynchrony, and increased interannual yield stability (Table 1). We found that the planted species richness treatments had prominent legacy effects on annually measured observed species richness throughout the 7-year duration of this study (Table 1) and that the negative relationship between observed species richness and weed abundance was an important driver of yield stability (Figs. 2 and 3). The following describes how our planted species richness and nitrogen fertility treatments affected the plant community through time—with an emphasis on relationships between planted species and weeds—and how those effects in turn influenced yield, yield variability, and, ultimately, interannual yield stability.

Table 1.

Analysis of variance results of response variables related to stability as affected by planted species richness, nitrogen, year, and their interactions.

Numbers in bold font indicate statistical significance (P < 0.05). P values followed by “+” indicate a positive and “−“indicate a negative relationship between main effects and response variables. Interactions are explained in the text.

	Yield*	Weed abundance†	Observed species richness*	Temporal variability of yield	Temporal stability‡	Species asynchrony
Planted richness (PR)	<0.001 +	<0.001 +	<0.001 +	0.605	0.002 +	<0.001 +
Nitrogen (N)	<0.001 +	0.004 −	<0.001 −	<0.001 +	0.861	0.715
Year (Y)§	<0.001 +	0.088	<0.001 +	NA	NA	NA
N × PR	0.182	0.244	0.026	0.987	0.273	0.984
PR × Y	<0.001	<0.001	0.049	NA	NA	NA
N × Y	0.148	0.019	0.044	NA	NA	NA
N × PR × Y	0.309	0.349	0.085	NA	NA	NA

*Response variable square root transformed.

†Results from analysis of deviance test for zero-inflated negative binomial model regression model.

‡Response log transformed.

§Year as a source of variability, and all interactions between year and other sources of variability, are not applicable (NA) for response variables estimated over time.

Fig. 2.

Structural equation model results.

Structural equation model results testing the relationships between observed species richness and nitrogen fertilization on species asynchrony and weed abundance, as well as their relationships between the two components of yield stability, yield and temporal variability of yield. Dashed lines indicate nonsignificant standardized path coefficients (P > 0.05), while red and blue solid lines indicate significant negative and positive (P < 0.05) standardized path coefficients, respectively. R2 values adjacent to each variable indicate the variance explained by the model for each respective variable.

Fig. 3.

Relationships between temporal stability and plant community characteristics.

Relationship between mean observed species richness (A), asynchrony (B), and weed abundance (C) on log-transformed temporal stability. Points are values from individual experimental units (n = 648) and colored coded by location (n = 9). The solid line indicates the estimated temporal stability value based on the x axis value as predicted by the linear mixed effects models described in Materials and Methods.

Species richness effects on yield, temporal variability of yield, and temporal stability

Increasing planted species richness directly reduced weed abundance [planted richness (PR): χ2 = 21.87; Table 1], and this effect was particularly dramatic in the high richness plots, where by the final year of the study, weed abundance was 96% less in 24 species mixtures than monocultures. Planted species richness was positively correlated with observed species richness (PR: F1,4404 = 9797.2; Table 1), which was used to identify and track relationships between the plant community and yield stability through time. Observed species richness was negatively associated with weed abundance (standardized path coefficient of direct effect, −0.20; Fig. 2), which was strongly negatively associated with biomass yield (standardized path coefficient of direct effect, −0.51). These linked relationships—observed species richness reducing weed abundance leading to subsequent yield increases—likely explain our important finding that effects of planted species richness on temporal mean yield are more influential to yield stability than temporal variability of yield (Fig. 4), and thus the positive relationship between planted species richness on yield stability (PR: F1,592 = 9.76; Table 1).

Fig. 4.

Relationship between planted species richness and both the mean yield (μ) and temporal variability of yield (σd).

Open circles represent mean yield values, and filled circles represent temporal variability of yield (g m−2) from unfertilized (A) and fertilized (B) experimental units (n = 324 for each fertilization treatment). Points are colored by location and jittered along the x axis to improve visualization. Solid lines represent estimated mean yield, and dashed lines represents estimated variability of yield from linear mixed effects models.

While increased plant diversity has been previously shown to reduce weed abundance in natural and agricultural settings (, ), this is, to our knowledge, the first study to link that relationship to temporal yield stability. We found that increasing species diversity reduced the negative impacts of weeds on mean biomass yield. Yield declines as a result of increased weed abundance may seem counterintuitive for grassland bioenergy systems because biomass from both native planted species and weed species contributes to yield. However, many weeds have short life cycles that do not correspond to those of native perennials and do not match harvest timings. Increasing diversity through the addition of weed species to existing perennial grassland plant communities does not result in increased productivity because weeds do not contribute to niche partitioning (i.e., decreased competition in diverse communities because intraspecific competition is stronger than interspecific competition), facilitation (i.e., increased facilitation in diverse communities because interspecific facilitation is stronger than intraspecific facilitation), and/or selection effects (i.e., the most productive species are also the most dominant species)—the common mechanisms used to explain how increased species diversity can increase productivity (). Our results emphasize the importance of limiting weed abundance to increase yields and temporal yield stability.

Another indirect pathway by which species richness was associated with reduced weed abundance and increased stability was through asynchrony (Figs. 2 and 3A). Species asynchrony is defined as the substitution, in terms of abundance, of species for others in a community in response to changing environmental conditions. Since plant communities often experience environmental variability through time, species asynchrony is hypothesized to be an important driver of the diversity-stability relationship (, ). We found a strong negative association between species asynchrony and weed abundance (standardized path coefficient of direct effect, −0.67), and its subsequent effects on yields were likely the main pathway for how diversity stabilized biomass yields (Figs. 2 and 3C). Communities with greater species richness and asynchrony can fill niche gaps that arise with interannual fluctuations in growing conditions, specifically niche gaps that are commonly filled by weed species. Filling niche gaps with planted native perennials can result in reduced weed abundance, and here, we show that this reduced weed abundance is correlated with higher mean yield.

Previous studies have shown that asynchrony can affect yield stability by reducing temporal variability of yield when species respond differently to year-to-year fluctuations in growing conditions (). Mixtures with asynchronous species that vary in functional traits [e.g., rooting depth () and N turnover ()] are more likely to maintain productivity in years with atypical growing conditions (). Our study did not include monocultures of all species represented in mixtures; thus, we were unable to confirm that compensatory responses to varying weather among species were directly responsible for the effect of asynchrony on yield stability. However, we did observe a positive correlation between functional group diversity and asynchrony (Pearson’s correlation coefficient, 0.31), indicating that plots with species varying by functional groups (C4 grasses, C3 grasses, legumes, and nonlegume forbs) did exhibit higher levels of asynchrony when experiencing changes in temperature and precipitation throughout the study (Fig. 1). Huang et al. () show that species and functional group asynchrony vary in their effect on plant community stability, especially in response to different types of environmental perturbations (e.g., warming and N enrichment). Selecting species for mixtures that vary in functional traits could help stabilize yields of bioenergy cropping systems.

Our study is the first to quantify relationships between weed abundance, planted species asynchrony, and yield stability. Our results are similar to () and () in that asynchrony was not positively associated with community productivity (standardized path coefficient of direct effect, −0.14), which we attribute to the dominant relationship between asynchrony and weed abundance, and subsequently weed abundance on mean yield (Fig. 2). Another notable finding here is that there was no significant relationship between asynchrony and temporal variability of yield—a relationship that has been previously observed and used to explain diversity–yield stability relationships (). Instead, species richness, asynchrony, and weed abundance influenced yield stability through a combination of direct and indirect effects on mean yield.

Nitrogen fertilization effects on yield, temporal variability of yield, and temporal stability

Nitrogen fertilization at 67 kg N ha−1 year−1 did not affect yield stability (N: F1,26 = 0.03) because it had opposite effects on the two components of stability: biomass yield and temporal variability of yield. Nitrogen fertilization increased yield (N: F1,26 = 84.56; Table 1; standardized path coefficient of direct effect, 0.34; Figs. 2 and 4) and increased temporal variability of yield (N: F1,26 = 30.28; Table 1; standardized path coefficient of direct effect, 0.24; Figs. 2 and 4). Averaged over planted species richness and years, fertilizer increased yields by 22%. This response to N fertilization was relatively low compared with the 52% increase in yield by native perennial mixtures and switchgrass monocultures when biomass yield was maximized at N fertilizer rates between 61 and 87 kg N ha−1 year−1 (). This low yield response to N averaged over species mixture treatments may be due to mixture-specific responses to N fertilizer. Although N fertilizer effects on yield did not vary by planted species richness in this analysis, Jungers et al. () reported that all treatments except the four-species legume mix and the eight-species grass/legume mix responded positively to N fertilization. A lack of yield response to N fertilization in these treatments may have diluted the average effect of N fertilization when averaged over all treatments. Unlike the effects of species richness, the effects of N fertilization were directly related to yield and temporal yield variability, and N fertilization was not associated with asynchrony and weakly related to weed abundance (standardized path coefficient of direct effect, 0.09).

Increased temporal variability in yield with N fertilizer could be related to the effects of N fertilization on species composition within the bioenergy mixtures and, subsequently, an increased dependence of such mixtures on other limiting resources [e.g., light; ()]. For example, annual fertilization of native plant communities shifts species composition so that species that are especially effective at acquiring and using N fertilizer become dominant (). This transition in species composition with N fertilization is often associated with an overall reduction in species richness (), which agrees with our finding that N fertilization decreased observed species richness (F1,26 = 18.34; Table 1). Since species that dominate under N fertilized conditions often do so by outcompeting neighboring species for light by growing taller (), these systems can result in high yields when conditions are especially suited for the dominant species. However, in years when precipitation or temperature limits the ability of the dominant species to reach maximum productivity, community yields can be especially low as there are fewer species in the mixture to compensate for the underperforming dominant species. Moreover, in years when environmental conditions limit assimilation of fertilizer N by N-dependent, dominant species—such as large spring precipitation events immediately following fertilization—detrimental environmental issues can occur including increased nitrate leaching to groundwater and/or nitrous oxide emissions (). In such cases, benefits of species diversity from both inter- and intra-annual niche partitioning can not only benefit yields but also limit environmental damage related to N fertilization.

Effects of diversity and species composition (identities of species) of plant communities can have strong effects on yield (), and fertilization can reduce diversity and shift composition in ways that reinforce () or counterbalance one another. The same is likely true for stability. Grman et al. () showed that N fertilization reduced species richness and resulted in species dominance, but, contrary to expectations, this transition improved yield stability. Their study was conducted at one site with coarse-structured, low organic matter soils and characterized as N limited. N fertilization eliminated the limiting resource, allowing the dominant species to thrive. Polley et al. () found that a single dominant species that was especially adapted to the range in interannual growing conditions was quite stable in terms of productivity, despite being dominant, and that this resulted in stable community yields. Ren et al. () also found that productivity by a dominant species was more important for community stability than diversity. The aforementioned studies (neither of which manipulated species richness or diversity as experimental treatments) show that dominant species can influence the relationships between diversity and stability, which highlights the value of long-term studies conducted at multiple sites like ours. Here, we have identified important relationships between plant community characteristics, N fertility, and productivity based on an experimental design that incorporates experimental units that may have had dominant species, were exposed to a broad range of temperature and precipitation patterns, and were tested across a range of soil texture classes and fertility levels (Fig. 1).

Plant community dynamics

Not all planted species established and/or persisted during our study. Although planted species richness was strongly correlated with observed species richness, the ratio of observed to planted species (slope parameter) was less in fertilized plots compared with unfertilized (N: F1,26 = 18.34; Table 1). Averaged over planted species richness and time, observed species richness was lower in fertilized (2.7 species) compared with unfertilized (2.9 species; N: F1,26 = 18.34; Table 1) treatments. Observed species richness slowly increased through time for unfertilized plots (0.06 species per year; N × Y: F1,4404 = 4.054; Table 1) but not for fertilized plots (0.02 species per year). Other studies have shown that the rate of species loss with N fertilization increases through time (), which will likely present negative indirect effects on yields that outweigh short-term positive effects. For short-term bioenergy production (<10 years), the negative effects of N fertilization on species richness would likely be accepted by bioenergy producers to achieve yield increases obtained with fertilization. However, producers may have to consider stand renovations or replacement for long-term maintenance of species richness and its positive effects on yield stability. Alternatively, producers could consider omitting N fertilization every few years or to reduce N fertilizer rates to maintain native species diversity.

The negative relationship between planted species richness and weed abundance increased through time (PR × Y: χ2 = 21.87; Table 1). Weed abundance in monocultures increased from 12.2 to 22.6% throughout the experiment, while weed abundance remained at 2.2% in the 24-species treatment averaged over time. Observed species richness also increased through time and more so in treatments with higher planted species richness (PR × Y: F1, 4404 = 3.86; Table 1). Higher observed species richness likely reduced above- and belowground resources available to weeds through time (). Observed species richness was negatively correlated with bare ground (Pearson’s correlation coefficient, −0.32), suggesting that space is one resource limited to weeds in species-rich mixtures. Removing one resource limitation—N through fertilization—did increase weed abundance, and this strengthened through time (N × Y: χ2 = 5.53; Table 1). At the onset of this study, N fertilization may have preferentially benefited weeds with superior N acquisition traits compared with the native species sown in the study. For example, annual weeds with higher relative growth rates compared with the perennial species may have become abundant early on, and reproduction and seed dispersal by these early colonizing weeds may have further increased weed abundance through time. Establishment of grasslands can be especially vulnerable to weed invasions when fertilized with N during the early years of establishment (), and increased weed abundance via N fertilization can increase the weed seedbank through time ().

Bioenergy yield potential (energy per unit area) is a function of biomass yield and feedstock composition, and plant species composition can influence both components of bioenergy potential. Biomass yield is the major source of variability in bioenergy potential (), and here, we show how plant species richness can stabilize biomass yield. Khalsa et al. () found that increasing species richness increased energy output of grassland systems harvested for combustion, yet some species and plant functional groups are known to contain especially high levels of elements problematic for gasification, such as silica in cool-season grasses (). Similarly, variability in fermentable carbohydrates important for cellulosic ethanol production has been documented by species and functional groups (). However, mixtures varying in native perennial species composition found in Minnesota did not vary in ethanol potential, and bioenergy yield potential was largely dependent on biomass yield ().

Observations and implications of diversity-productivity relationships

Under the conditions of this study—where species composition was not managed, weeds were judiciously controlled, and biomass was harvested annually—we observed a positive relationship between planted species richness and biomass yield. Few studies have shown a similar relationship under such conditions (, ), but these findings are supported by many grassland diversity experiments around the globe (, , ). However, here we report some differences with some other diversity experiments. First, the positive effects of species richness on yield weakened through time (PR × Y: F1,3972 = 12.78) in our study. This is contrary to findings from most grassland diversity experiments () including those ranging in durations from four () to more than 13 years (). This temporal trend observed in our study is likely related to relatively low yields of the monoculture grass species during early years of establishment, whereas the 12 and 24 species mixtures had similar yields early in the experiment compared with those at the end (). Therefore, the difference in yield between monocultures and diverse mixtures was greater in the early years of the study but became more similar through time. For producers planning to manage stands for less than 10 years, this finding may have little consequence on production and profitability. However, producers considering long-term management of perennial bioenergy crops will need to consider methods other than plant diversity to maintain yields as stands age.

The economic outcome of increasing species diversity in perennial grassland bioenergy mixtures is dependent on highly variable and uncertain prices of inputs (seed) and outputs (biomass). Some of the species included in the 24–species mixture treatment are relatively uncommon and have not been considered for bioenergy production because of their low biomass potential and relatively high seed cost. Economic models that can simulate various scenarios of highly variable seed and biomass prices, as well as incorporate potential indirect benefits related to plant diversity (e.g., increased soil organic matter) and direct benefits in the form of payments for ecosystem services, could help producers make economically sound decisions on species selection for bioenergy crop production. Moreover, future experiments conducted using commercial-scale plots and equipment can provide more accurate estimates of real-world yield and economic outcomes ().

The detrimental effects of weeds on annual food crop production are well known. However, effects of weeds on perennial bioenergy crops are less obvious since weed biomass is perceived to function and contribute similarly to bioenergy yields as any other species in a polyculture mixture. We quantified the negative relationship between weed abundance on biomass yields and relate that to a reduction in yield stability. We also show that increasing native species richness in perennial polyculture bioenergy cropping systems promotes biomass yield and stability by reducing weed abundance. These findings illustrate new ways that increasing diversity can increase the economic potential of a perennial bioenergy cropping systems through interannual stability in yields. Yield stability and the economic security it provides can be expected to become increasingly important as precipitation and temperature become more variable with anticipated climate change.

MATERIALS AND METHODS

Experimental design

The experiment was conducted at nine sites located in Minnesota and North Dakota spanning a large gradient of precipitation, temperature, and soil types (Fig. 1; Supplementary Materials, table S1). Details on climate history and soil types can be found in () and (). The crop productivity index (CPI) was estimated for each site (Fig. 1) based on the National Commodity Crop Productivity Index model developed by the U.S. Department of Agriculture National Resource Conservation Service. The CPI ratings range from 0 to 100, with higher ratings indicating greater maize (Zea mays L.) yield potential and are based on physical and chemical soil properties from the National Cooperative Soil Survey data. Sites north of 46°N were seeded in 2007 and south of 46°N were seeded in 2006. The experiment was a randomized, split-plot design with three replications per location. The whole-plot treatment was N fertilizer rate applied at two levels, either 0 or 67 kg N ha−1 year−1 broadcast in spring in the form of ammonium nitrate. Split-plots (each 9 m2 in size) were randomly assigned to one of 12 species mixture treatments with varying levels of plant diversity from 1, 4, 8, 12, or 24 species that were purposefully selected for productivity and regional suitability (SI Appendix and table S2). The species selected consisted of C4 grasses, C3 grasses, legumes, and nonleguminous forbs. All treatments were seeded at a rate of 12 g m−2. For mixture treatments, seed mass was divided equally among the species assigned to the treatments. A complete list of species and seeding rates can be found in (). Nitrogen fertilizer treatments were first imposed in 2008, one or two growing seasons after sowing, and then applied annually for all remaining years of the study. Weeds were managed using tools and practices available for commercial production of biomass feedstock from perennial polyculture grasslands composed of monocot and dicot species. Before 2008, minimal weed control was conducted so that experimental plots resembled species composition and yield potential of fields under large-scale production. Herbicides were not used in the first 2 years to avoid injury to the establishing seedlings that emerged from sown seeds. Only rarely were weeds removed by hand when the weed canopy closed and excluded light from sown seedlings. Acetochlor [2-chloro-N-(ethoxymethyl)-N-(2-ethyl-6-methylphenyl) acetamide] was applied at all locations in late April in 2008 and 2009, but application was forgone after poor weed control was observed. In subsequent years, plots were largely unweeded, except in cases of local outbreaks of Canada thistle [Cirsium arvense (L.) Scop.], which occurred at some sites and were controlled by spot spraying of glyphosate and/or removal by hand.

The relative abundance of species was assessed by visually estimating the percent ground cover of each planted species shortly after peak biomass. During this survey, the relative abundance of weeds was also estimated as percent ground cover. Weeds were not identified to the species level and did not include volunteer native species that may have colonized from neighboring plots (which was minimal). Exposed soil was categorized as bare ground. The same observer measured species ground cover at all locations and across all years of the study. Biomass yield was harvested annually following a killing frost (−2°C). Within each subplot, all biomass was clipped manually to a stubble height of ~5 cm from a 0.75-m2 sample quadrat. Harvested biomass was then dried at 60°C until constant mass before weighing for dry matter yield.

There were 72 experimental units per location (12–species mixture treatments × 2 N fertilizer rates × 3 replications). With 7 years of data from four southern locations (2007 to 2013) and 6 years of data from five northern locations (2008 to 2013), our dataset included 4176 total observations of biomass yield and percent cover of planted species and weeds. Temporal stability (μ/σd) was calculated as the mean yield divided by the detrended temporal variability of yield for each experimental unit (). Mean yield (μ) was calculated as the average of biomass yields across years for each experimental unit. Temporal variability was calculated following () to account for general trends in yield through time. Instead of using the SD of biomass yield through time for each experimental unit, a linear model with yield as the response variable and year as a predictor was constructed for each experimental unit, and the residual standard error of the model was used as a detrended statistic for temporal variability (σd). Species synchrony was calculated for each experimental unit following () using the synchrony function in the “vegan” R package (). Asynchrony was calculated by subtracting the synchrony value from 1. Species asynchrony was calculatedwith Y as the abundance of species i in a plot of n species, and Y is the abundance of all other species in the community except Y. Functional group diversity was calculated using the Shannon diversity index and the relative proportion of C4 grasses, C3 grasses, legumes, and nonlegume forbs from species percent cover assessments for all plots assigned to a species mixture treatment.

Statistical analysis

Linear mixed effects models were used to partition variation in annual biomass yield, weed abundance, and observed species richness into sources including planted species richness (continuous), N fertilizer (categorical), year since establishment (continuous), and all possible two- and three-way interactions (). Random effects included N fertilizer treatment nested within block nested within location to account for pseudoreplication of split plot treatments. To meet model assumptions, yield and observed species richness data were square root transformed, and weed abundance data were modeled using zero-inflated negative binomial regression. Linear mixed effects models were used to partition variation in temporal variability of biomass yield, temporal stability, and species asynchrony into sources including planted species richness (continuous), N fertilizer (categorical), and all their two-way interactions. The same random effects structure was used as described above. Temporal stability was log transformed to meet model assumptions. All analyses were conducted with R version 3.5.3.

We used a structural equation modeling approach to understand the direct and indirect effects of species richness and nitrogen fertilizer application on stability. We specified the models using species richness, N application, asynchrony, weed cover, temporal variability of yield, and mean yield as parameters explaining stability. Data were fitted using the maximum likelihood estimation method with the package “lavaan” in R. Model fit was determined using χ2, root square mean errors of approximation, and Akaike information criterion. Last, standardized parameter estimates were used to compare the effect sizes of the parameters in the final model.