Letter doi:10.

1038/nature16524

Integrative modelling reveals mechanisms linking

productivity and plant species richness
James B. Grace1, T. Michael Anderson2, Eric W. Seabloom3, Elizabeth T. Borer3, Peter B. Adler4, W. Stanley Harpole5,6,7,
Yann Hautier8, Helmut Hillebrand9, Eric M. Lind3, Meelis Pärtel10, Jonathan D. Bakker11, Yvonne M. Buckley12,
Michael J. Crawley13, Ellen I. Damschen14, Kendi F. Davies15, Philip A. Fay16, Jennifer Firn17, Daniel S. Gruner18,
Andy Hector19, Johannes M. H. Knops20, Andrew S. MacDougall21, Brett A. Melbourne15, John W. Morgan22, John L. Orrock14,
Suzanne M. Prober23 & Melinda D. Smith24

How ecosystem productivity and species richness are interrelated is The search for a canonical bivariate productivity–richness relation-
one of the most debated subjects in the history of ecology1. Decades ship lies at the heart of the debate among ecologists. This pursuit is
of intensive study have yet to discern the actual mechanisms behind fuelled, in part, by the history of the discussion, which has focused
observed global patterns2,3. Here, by integrating the predictions on bivariate predictions16. At the same time, it is also seen by some as
from multiple theories into a single model and using data from a means of assessing the overall importance of various mechanisms
1,126 grassland plots spanning five continents, we detect the clear operating in natural systems. While many different mechanisms have
signals of numerous underlying mechanisms linking productivity been discussed, the primary competing theories make four main
and richness. We find that an integrative model has substantially conflicting predictions: (1) richness and productivity should increase
higher explanatory power than traditional bivariate analyses. In together with increasing resources and environmental favourabil-
addition, the specific results unveil several surprising findings that ity until limits to coexistence are reached at high productivity and
conflict with classical models4–7. These include the isolation of a richness declines, producing a humped-shape relationship4–7,17;
strong and consistent enhancement of productivity by richness, (2) richness promotes productivity, leading to a positive relationship6,9;
an effect in striking contrast with superficial data patterns. Also (3) richness and productivity increase together because climatic gradi-
revealed is a consistent importance of competition across the ents in productivity lead to increased regional species pools, creating
full range of productivity values, in direct conflict with some a positive relationship but from a separate mechanism8; and (4) the
(but not all) proposed models. The promotion of local richness richness–productivity relationship will be of inconsistent form because
by macroecological gradients in climatic favourability, generally the mechanisms controlling them vary in their scale-dependence and
seen as a competing hypothesis8, is also found to be important in relative importance18,19.
our analysis. The results demonstrate that an integrative modelling Empirical tests of the generality of hypothesized bivariate produc-
approach leads to a major advance in our ability to discern the tivity–richness patterns have reported a wide variety of results and
underlying processes operating in ecological systems. have produced substantial discussion20–23. Recent global studies2,3
Ecosystem productivity and species diversity are essential to the have disagreed with regard to whether a coherent pattern exists for
ability of natural systems to provide goods and services. Yet, for dec- natural grasslands. What has been agreed upon, however, is that the
ades there has been debate over their interrelationship. In the 1970s low explanatory power coming from conventional analyses suggests
and 1980s, conflicting models predicted that elevated productivity the need to pursue an integrative understanding of the causal mecha-
would lead to reductions in species richness4–7. Beginning in the mid- nisms controlling productivity–richness relationships.
1990s, scientists started to seriously debate another possibility: that One potential explanation for why debate over mechanisms is prov-
richness could promote productivity9–12. While experimental studies ing difficult to resolve is because productivity and richness are jointly
generally support a biodiversity enhancement of productivity13,14, controlled by a complex network of processes1,21,24–28. Overcoming
the precise strength of the effect in natural systems and the relation- the challenge of evaluating more complex hypotheses requires both
ship of this process to other factors that can influence productiv- advanced statistical modelling approaches and large-scale systematic
ity remain major questions. Adding to the debate, macroecological data collection efforts. Here, we used structural equation modeling29 to
theories propose that regional diversity is controlled by gradients integrate key predictions from competing theories into a multi-process
in climatic favourability and evolutionary history15 and that these hypothesis for evaluation. We then evaluated the hypothesis using
larger-scale effects are important determinants of smaller-scale diver- data collected for that purpose by a global consortium, the Nutrient
sity patterns8. Network ([Link] The data collected comprise samples
1
US Geological Survey, Wetland and Aquatic Research Center, 700 Cajundome Boulevard, Lafayette, Louisiana 70506, USA. 2Department of Biology, 206 Winston Hall, Wake Forest University,
Box 7325 Reynolda Station, Winston-Salem, North Carolina 27109, USA. 3Ecology, Evolution, and Behavior, University of Minnesota, 1987 Upper Buford Circle, St Paul, Minnesota 55108,
USA. 4Department of Wildland Resources and the Ecology Center, Utah State University, 5230 Old Main, Logan, Utah 84322, USA. 5Department of Physiological Diversity, Helmholtz Center for
Environmental Research – UFZ, Permoserstrasse 15, 04318 Leipzig, Germany. 6German Centre for Integrative Biodiversity Research (iDiv), Deutscher Platz 5e, D-04103 Leipzig, Germany.
7
Martin Luther University Halle-Wittenberg, Am Kirchtor 1, 06108 Halle (Saale), Germany. 8Ecology and Biodiversity Group, Department of Biology, Utrecht University, Padualaan 8, Utrecht
3584 CH, The Netherlands. 9Institute for Chemistry and Biology of the Marine Environment, University of Oldenburg, Schleusenstrasse 1, Wilhelmshaven D-26381, Germany. 10Institute of Ecology
and Earth Sciences, University of Tartu, Lai 40, Tartu 51005, Estonia. 11School of Environmental and Forest Sciences, University of Washington, Box 354115, Seattle, Washington 98195-4115,
USA. 12School of Natural Sciences, Zoology, Trinity College Dublin, The University of Dublin, Dublin 2, Ireland. 13Department of Biological Sciences, Imperial College London, Silwood Park, Ascot,
Berkshire SL5 7PY, UK. 14Department of Zoology, University of Wisconsin, 430 Lincoln Drive, Madison, Wisconsin 53706, USA. 15Department of Ecology and Evolutionary Biology, UCB 334,
University of Colorado, Boulder, Colorado 80309, USA. 16Grassland Soil and Water Research Laboratory, United States Department of Agriculture Agricultural Research Service, 808 East Blackland
Road, Temple, Texas 76502, USA. 17#15 Queensland University of Technology, School of Earth, Environment and Biological Sciences, Brisbane, Queensland 4001, Australia. 18Department of
Entomology, University of Maryland, College Park, 4112 Plant Sciences, College Park, Maryland 20742, USA. 19Department of Plant Sciences, University of Oxford, South Parks Road, Oxford
OX1 3RB, UK. 20School of Biological Sciences, 348 Manter Hall, University of Nebraska, Lincoln, Nebraska 68588, USA. 21Department of Integrative Biology, University of Guelph, Guelph, Ontario
N1G 2W1, Canada. 22Department of Ecology, Environment, and Evolution, La Trobe University, Bundoora, Victoria 3083, Australia. 23CSIRO Land and Water, Private Bag 5, Wembley, Western
Australia, 6913, Australia. 24Department of Biology, Colorado State University, 1878 Campus Delivery, Fort Collins, Colorado 80526, USA.

0 0 M o n t h 2 0 1 6 | VO L 0 0 0 | NAT U R E | 1
© 2016 Macmillan Publishers Limited. All rights reserved
 RESEARCH Letter

A Bivariate plot Figure 1 | Comparison between low-dimension

60
(top panel) and high-dimension (bottom panel)
400
examinations of data. A, Raw bivariate plot of

Frequency
above-ground productivity and species richness
50
in 1-m2 plots (n = 1,126). Different sites are
200
Species richness

0 represented in the graphs by different colours,

40 0 500 1,000 1,500 assigned by mean site richness from low (yellow)
Plot productivity
to high (red). B, Plots a–c visualize site level partial
30 400
relationships indicated by corresponding letters in
Fig. 2 (n = 39 sites). Plots d–f visualize plot level

Frequency
20 200 partial relationships indicated by corresponding
letters in Fig. 2 (n = 1,126 1-m2 plots). Units are
0 0 standardized residual deviations from predicted
0 10 20 30 40 partial scores.
0 500 1,000 1,500 Plot richness

Annual productivity (g)

B Multivariate partial plots

a b c
2
2
1
Site productivity
Site richness

Site richness
1
1
0
0
0
−1 −1
−1
−2 −2
−2

−1 0 1 2 3 −1 0 1 2 3 −1.5 −0.5 0.5 1.5

Site biomass Site richness Mean annual precipitation
d e f 4
2 2
1
Plot productivity

1 2
Plot richness
Plot shading

0
0 0
−1
−1 −2 −2
−2 −3
−4
−3 −2 −1 0 1 2 −1.5 −0.5 0.5 1.5 −4 −2 0 2
Plot biomass Plot shading Plot richness

from 1,126 plots collected at 39 grass-dominated sites around the model reveals strong, clear signals consistent with numerous proposed
world. Variables measured include plant species richness, productivity mechanisms, including several that are not at all suggested from the
(measured as the annual biomass increment), total biomass (the accu- bivariate data.
mulated non-woody biomass, live and dead, including litter), along First, we found clear evidence that the accumulation of total bio-
with many of the drivers hypothesized to be important for regulating mass (hereafter simply ‘biomass’) leads to a negative effect on spe-
their variations. Additional information is provided in the Methods. cies richness. At the site level, the partial effect (r∂) of biomass on
To integrate theoretical expectations from competing theories, we richness in the model was strong (Figs 1B, a and 2; r∂ = −0.77). The
mined the productivity–diversity literature to determine the main reduction of richness was not found to be mediated by our one-
theoretical constructs discussed and the hypothesized interconnec- time measurement of average shading at the ground surface, which
tions between constructs (see Methods). We used this information to was subsequently dropped from the model. At the plot level, how-
develop a structural equation meta-model that assimilates the essential ever, we found evidence that biomass increases shading (r∂ = 0.56),
theoretical constructs and hypothesized connections into a network of which in turn, decreases richness (Figs 1B, d, e and 2; r∂ = −0.34).
multivariate expectations (Extended Data Fig. 1, Extended Data Table 1 The negative effects of biomass on richness appear consistent with
and Supplementary Information). This meta-model, along with the long-standing hypotheses that predict a hump-shaped productivity–
available data, guided our development of a structural equation model richness relationship due to competitive dominance at high pro-
for empirical evaluation. We evaluated model-data consistency to deter- ductivity5,17. However, while those hypotheses assume increasing
mine whether there were missing linkages in the initial model as well as competitive intensity with increasing productivity, our results reveal
to determine the support for proposed links. We further addressed the a linear effect across the full range of biomass observed in this study
question, ‘what dimension of model (that is, number of parameters and (Fig. 1B, a).
linkages) is required to detect the signals in the underlying data?’. For Second, we found a positive, linear enhancement of productivity
this, we evaluated lower-dimensional versions of the model by remov- by richness in the model. This effect was among the strongest found
ing linkages and re-evaluating against the data. More methodological at the site-scale (Figs 1B, b and 2; r∂ = 0.67), and was detectable,
detail is provided in the Methods and Supplementary Information. although weak, at the plot-scale (Figs 1B, f and 2; r∂ = 0.02). A surpris-
A simple bivariate plot of richness against productivity (Fig. 1A) ing feature of the site-level result is the apparent absence of a levelling
reveals little about the underlying mechanisms. Previous analyses of off of the biodiversity enhancement of productivity at higher levels
such bivariate relations have found it difficult to even detect signif- of richness. Such a continuous effect has been theorized for larger-
icant associations2,30. However, our analysis based on an integrative scale studies and contrasts with the asymptotic levelling off usually

2 | NAT U R E | VO L 0 0 0 | 0 0 m o n t h 2 0 1 6
© 2016 Macmillan Publishers Limited. All rights reserved
 Letter RESEARCH

Table 1 | Select standardized partial effect sizes (and standard errors) ranked by magnitude and proposed interpretations
Effect Magnitude Proposed interpretation

Site-scale
Soil fertility → productivity 1.104 (0.220) Productivity variations are most strongly related to spatial variations in site fertility.
Site biomass → richness −0.771 (0.143) High biomass sites have depressed species richness, presumably via some form of competitive dominance.
Site richness → productivity 0.671 (0.214) Increasing richness contributes to higher productivity.
Climate → richness 0.669 (0.113) Richness increases with increasing mean annual precipitation during warmest quarter.
Heterogeneity → richness 0.627 (0.119) Species coexistence strongly regulated by within-site spatial heterogeneity in vegetation cover.
Climate → productivity 0.589 (0.195) Combined effects of macrogradients in temperature and precipitation are moderately important in
controlling site-to-site variations in productivity.
Soil suitability → richness 0.439 (0.107) Richness increases with decreasing sand and silt content of soils.
Disturbance → richness −0.251 (0.116) Sites subject to strong anthropogenic alteration (for example, pasturing) notably lower in richness.
Disturbance → biomass −0.185 (0.084) Local grazing contributed to some biomass variation.
Correlation between soil fertility −0.56 (0.109) Influences of resources and filters on richness are distinct from those that regulate productivity.
and soil suitability
Plot-scale
Biomass → shading 0.559 (0.089) Shading fairly strongly tied to biomass, although morphology and other features probably also important.
Soil suitability → richness 0.404 (0.074) Plot-to-plot variations in richness tied to local variations in environmental edaphic filters.
Shading → richness −0.342 (0.064) Richness variations within site partially controlled by variations in competition via shading.
Soil suitability → shading 0.249 (0.098) Variations in community traits driven by soil resources result in differences in amount of shading per gram.
Richness → productivity 0.017 (0.006) Effects of natural plot-to-plot variations in richness can be discerned as an independent process in this
sample.
Soil fertility → productivity 0.0 (NS) Plot productivity largely determined by site-level productivity, with no significant local control.
Site-scale and plot-scale effects are presented separately.

found in experimental (smaller-scale) studies13. Previous attempts correlated at the site level (Table 1, r∂ = −0.56), supporting previous
to isolate an effect of richness on productivity with observational claims of their semi-independence21. Thus, theories that presume a
data using simpler models have failed to do so (see Supplementary simultaneous increase in productivity and richness with increasing
Information). environmental favourability (see Supplementary Information) fail to
Third, we found strong and independent influences of macroclimate correspond with the independent responses to environmental drivers
and soils on richness and productivity. The standardized effect sizes pro- observed in natural systems.
vide insights into the relative importance of these processes (Table 1). Our results show that failure to account for the variation in richness
At the site level, productivity was most strongly related to soil fertil- and productivity explained by the environmental drivers would make
ity, while richness was most strongly related to climate and soil suit- it difficult to detect the reciprocal influences of productivity and rich-
ability, with heterogeneity and disturbance also important (Fig. 2). ness on each other. In fact, our capacity to isolate underlying processes
Rather than being made up of similar environmental factors, the soil was highly sensitive to model dimensionality, where dimensionality
environmental drivers of richness and productivity were negatively refers to the number of measured determinants of productivity and

Climatic
Figure 2 | Structural equation model
resources, Disturbance Heterogeneity
representing connections between productivity
regulators, and richness supported by the data. ‘Biomass’
and filters
0.67 –0.25
refers to total above-ground accumulated
–0.18 biomass. Letters correspond to partial plots
0.59 shown in Fig. 1B. Solid arrows represent positive
0.63
effects, dashed arrows represent negative effects.
c For the site-level submodel, test statistic = 13.518,
with 13 model degrees of freedom and P = 0.409
Site 1.1 Site –0.77 a Site (indicating close model-data fit). For the plot-
productivity
2
biomass
2
richness level submodel, robust test statistic = 21.907, with
R = 0.36 R = 0.73 R2 = 0.61
16 model degrees of freedom and P = 0.146 (again
indicating close model-data fit). Relative effect
sizes presented in Table 1.
0.67 b
0.85 0.67

0.89 0.02 f

1.1
0.44
Plot Plot d e Plot
Shading
productivity biomass richness
2 NS 2 0.56 R2 = 0.43 2
R = 0.72 R = 0.79 –0.34 R = 0.65

NS 0.40

Soil Soil
fertility suitability

0 0 M o n t h 2 0 1 6 | VO L 0 0 0 | NAT U R E | 3
© 2016 Macmillan Publishers Limited. All rights reserved
 RESEARCH Letter

richness included in the model (Extended Data Table 2). At both site 8. Gillman, L. N. & Wright, S. D. The influence of productivity on the
and plot levels, models omitting either productivity or biomass (but species richness of plants: a critical assessment. Ecology 87, 1234–1243
(2006).
not both) still permitted us to detect the feedback from richness to 9. Naeem, S., Thompson, L. J., Lawler, S. P., Lawton, J. H. & Woodfin, R. M.
biomass production. Any other simplifications at the site-level, how- Declining biodiversity can alter the performance of ecosystems. Nature
ever, resulted in a failure to detect previously detected pathways and 368, 734–737 (1994).
10. Tilman, D., Wedin, D. & Knops, J. Productivity and sustainability influenced by
resulted in a dramatic loss of signal (as indicated by reduced values biodiversity in grassland ecosystems. Nature 379, 718–720 (1996).
of R2 in the model). 11. Huston, M. A. Hidden treatments in ecological experiments: re-evaluating the
Regarding scale dependence, plot-level values of productivity, bio- ecosystem function of biodiversity. Oecologia 110, 449–460 (1997).
12. Grime, J. P. Biodiversity and ecosystem function: the debate deepens. Science
mass, and richness were strongly related to site-level estimates (Fig. 2), 277, 1260–1261 (1997).
as is common with hierarchical data. This should be interpreted as 13. Cardinale, B. J. et al. Biodiversity loss and its impact on humanity. Nature 486,
meaning much of the overall plot-to-plot variation in productivity, 59–67 (2012).
14. Tilman, D., Reich, P. B. & Isbell, F. Biodiversity impacts ecosystem productivity
biomass, and richness can be ascribed to site-to-site variations in those as much as resources, disturbance, or herbivory. Proc. Natl Acad. Sci. USA 109,
properties. In this case, within-site variations in productivity were 10394–10397 (2012).
explained solely by site-level productivity, as there were no predictors 15. Hawkins, B. A. et al. Energy, water, and broad-scale geographic patterns of
for remaining among-plot variations. Within-site variations in rich- species richness. Ecology 84, 3105–3117 (2003).
16. Pärtel, M., Zobel, K., Laanisto, L., Szava-Kovats, R. & Zobel, M. The productivity-
ness, however, were additionally explained by within-site variations in diversity relationship: varying aims and approaches. Ecology 91, 2565–2567
soil suitability and shading. Also sensitive to scale was the strength of (2010).
the feedback from richness to productivity, which was much stronger 17. Grime, J. P. Plant Strategies and Vegetation Processes (John Wiley, 1979).
18. Chase, J. M. & Leibold, M. A. Spatial scale dictates the productivity-biodiversity
at the site scale. While multiple factors probably play a role in this relationship. Nature 416, 427–430 (2002).
scale-dependence, the simplest explanation here may be the smaller 19. Scheiner, S. M. & Jones, S. Diversity, productivity and scale in Wisconsin
span of conditions sampled within sites compared with across sites. vegetation. Evol. Ecol. Res. 4, 1097–1117 (2002).
20. Waide, R. et al. The relationship between productivity and species richness.
Finally, in contrast to a bivariate model, which our analyses suggest Annu. Rev. Ecol. Syst. 30, 257–300 (1999).
can explain no more than 10% of the observed variation in richness, 21. Grace, J. B. The factors controlling species density in herbaceous plant
our structural equation model explains 61% of the variation in richness communities: an assessment. Perspect. Plant Ecol. Evol. Syst. 2, 1–28
(1999).
among sites, and 65% of the variation in richness among plots. An 22. Mittelbach, G. G. et al. What is the observed relationship between species
ability to explain a substantial portion of the variation in richness is richness and productivity? Ecology 82, 2381–2396 (2001).
tremendously important for potential conservation applications. Model 23. Whittaker, R. J. Meta-analyses and mega-mistakes: calling time on
meta-analysis of the species richness-productivity relationship. Ecology 91,
complexity is also important because of its more detailed mapping onto 2522–2533 (2010).
nature, as our model can make statements about how both specific 24. Borer, E. T. et al. Herbivores and nutrients control grassland plant diversity via
management actions (such as reduction of biomass through mowing light limitation. Nature 508, 517–520 (2014).
or increase in soil fertility through fertilization), as well as shifts in 25. Schmid, B. The species richness–productivity controversy. Trends Ecol. Evol. 17,
113–114 (2002).
climate conditions, may alter both productivity and species richness. 26. Loreau, M., Naeem, S. & Inchausti, P. Biodiversity and Ecosystem Functioning:
Our findings give reason for optimism about the future of ecology Synthesis and Perspectives (Oxford Univ. Press, 2002).
as a more precise and less ambiguous science. We show that many of 27. Cardinale, B. J., Bennett, D. M., Nelson, C. E. & Gross, K. Does productivity drive
diversity or vice versa? A test of the multivariate productivity-diversity
the proposed processes connecting productivity and richness offered hypothesis in streams. Ecology 90, 1227–1241 (2009).
during previous decades operate simultaneously as parts of a whole 28. Scherber, C. et al. Bottom-up effects of plant diversity on multitrophic
system of effects. Details of the findings, however, refine many of our interactions in a biodiversity experiment. Nature 468, 553–556 (2010).
29. Grace, J. B. et al. Guidelines for a graph-theoretic implementation of structural
assumptions about how those processes operate. Our field’s previous equation modeling. Ecosphere 3, art73 (2012).
failure to resolve debate about productivity–richness relationships 30. Grace, J. B., Adler, P. B., Stanley Harpole, W., Borer, E. T. & Seabloom, E. W.
stems from a lack of integration of ideas and absence of simultane- Causal networks clarify productivity–richness interrelations, bivariate plots do
not. Funct. Ecol. 28, 787–798 (2014).
ous tests of their combined implications. By integrating and testing
those ideas, our approach provides a systems-level understanding and Supplementary Information is available in the online version of the paper.
improves our chances to foresee the possible consequences of human
Acknowledgements J.B.G. was supported by the US Geological Survey
alteration of environmental factors, productivity, and richness now Ecosystems and Climate and Land use Change Programs. This work uses
occurring worldwide. data from the Nutrient Network ([Link] experiment, funded at
the site scale by individual researchers. Coordination and data management
Online Content Methods, along with any additional Extended Data display items and were supported by funding to E.T.B. and E.W.S. from the National Science
Source Data, are available in the online version of the paper; references unique to Foundation (NSF) Research Coordination Network (NSF-DEB-1042132) and
these sections appear only in the online paper. Long Term Ecological Research (NSF-DEB-1234162 to Cedar Creek LTER)
programs and the UMN Institute on the Environment (DG-0001-13). The
received 1 September; accepted 8 December 2015. Minnesota Supercomputer Institute hosts project data. The use of trade,
firm, or product names is for descriptive purposes only and does not imply
Published online 13 January 2016. endorsement by the US Government. Support for site-level activities is
acknowledged in the Supplementary Information. We thank D. Laughlin for
1. Willig, M. R. Biodiversity and productivity. Science 333, 1709–1710 (2011). comments on the manuscript.
2. Adler, P. B. et al. Productivity is a poor predictor of plant species richness.
Science 333, 1750–1753 (2011). Author Contributions E.W.S., E.T.B., W.S.H. and E.M.L. are Nutrient Network
3. Fraser, L. H. et al. Plant ecology. Worldwide evidence of a unimodal relationship coordinators. J.B.G. and T.M.A. developed and framed the research questions.
between productivity and plant species richness. Science 349, 302–305 T.M.A., E.W.S., E.T.B., P.B.A., W.S.H., Y.H., H.H., J.D.B., Y.M.B., M.J.C., E.I.D., K.F.D.,
(2015). P.A.F., J.F., D.S.G., A.H., J.M.H.K., A.S.M., B.A.M., J.W.M., J.L.O., S.M.P. and M.D.S.
4. Grime, J. P. Competitive exclusion in herbaceous vegetation. Nature 242, collected data used in this analysis. T.M.A. assembled the data and performed
344–347 (1973). initial analyses. J.B.G. analysed the data and wrote the paper with contributions
5. Huston, M. A. A general hypothesis of species diversity. Am. Nat. 113, 81–101 and input from all authors.
(1979).
6. Tilman, D. Resource Competition and Community Structure (Princeton Univ. Author Information Reprints and permissions information is available at
Press, 1982). [Link]/reprints. The authors declare no competing financial
7. Taylor, D. R., Aarssen, L. W. & Loehle, C. On the relationship between r/K interests. Readers are welcome to comment on the online version of the paper.
selection and environmental carrying capacity: a new habitat templet for plant Correspondence and requests for materials should be addressed to
life history strategies. Oikos 58, 239–250 (1990). J.B.G. (gracej@[Link]).

4 | NAT U R E | VO L 0 0 0 | 0 0 m o n t h 2 0 1 6
© 2016 Macmillan Publishers Limited. All rights reserved
 Letter RESEARCH

Methods Laboratory, Memphis, Tennessee, USA; these included the following: extractable
Development of meta-model hypothesis. A review and accounting of the history soil phosphorus and potassium were quantified using the Mehlich-3 extraction
of claims and disputed points in the published literature was developed before method, and parts per million concentration estimated using inductively coupled
construction of the meta-model that guided this analysis (Extended Data Fig. 1 plasma-emission spectrometry. Soil pH was quantified with a pH probe (Fisher
and Supplementary Information). During this review, attention was paid to the Scientific) in a slurry made from 10 g dry soil and 25 ml of deionized water. Soil
theoretical constructs invoked by various authors, since our goal was to provide a texture, expressed as the percentage sand, percentage silt, and percentage clay,
framework that had the potential to clarify and resolve disputed points. Attention was measured on 100 g dry soil using the Buoycous method. Further details on
was also paid to types of variable measured by different authors, as the relationship sampling methodology are at [Link]
between constructs and measurements constitutes one of the several sources of Climatic characteristics were obtained for each site from version 1.4 of BioClim,
ambiguity and confusion31,32. An in-depth description of the literature synthesized which is part of the WorldClim40 set of global climate layers at 1 km2 spatial reso-
to generate the meta-model is presented in the Supplementary Information. lution. To represent measures of temperature and precipitation with meaningful
Data. Data collected by the Nutrient Network Cooperative33 was used to design relationships to plant growth in global grasslands, we selected mean temperature of
and evaluate a structural equation model based on the meta-model presented. The the wettest quarter of the year (BIO8) and total precipitation of the warmest quar-
Nutrient Network is a distributed, coordinated research cooperative. Sites in the ter of the year (BIO18). Climate values were extracted using universal transverse
Network are dominated primarily by herbaceous vegetation and intended to repre- Mercator (UTM) coordinates collected near the centre of each site.
sent natural/semi-natural grasslands and related ecosystems worldwide. Individual Several derived variables were developed to include in the modelling effort. To
sites were selected to accommodate at least a 1,000 m2 study design footprint. Most represent within-site heterogeneity, coefficients of variation were computed for the
sites sampled vegetation in 2007, although 12 sites sampled in 2008 or 2009. No site-level model based on plot-to-plot variation in plot-level measures. This allowed
statistical methods were used to predetermine sample size. Samples were collected us to examine the explanatory value of heterogeneity in soil nitrogen, phosphorus,
using a completely randomized block design. The standard design has three blocks potassium, and pH, as well as heterogeneity in biomass and light interception.
and ten plots per block at each site, although some sites deviate slightly from this Indices of total resource supply and resource imbalance were also calculated using
design. A few sites are grazed or burned before sampling, consistent with their the method of ref. 27 and evaluated for inclusion in our models.
traditional management. Further details on site selection and design can be found Disturbance history information for the sites was converted into four binary
at [Link] (0,1) variables for analyses; information available included pretreatment history
In this study, we analysed data from 39 of the 45 sites considered in ref. 2 pos- of (1) substantial anthropogenic alteration (for example, conversion to pasture),
sessing a complete set of covariates (Extended Data Table 3). While ref. 2 only (2) grazing history, by wild or domestic animals, (3) active management (typically
examined bivariate relations between productivity and richness, our analyses haying or mowing), and (4) fire. Current levels of herbivory were estimated by
brought in many additional variables (Extended Data Table 1) so that we could comparing biomass inside and outside exclosure plots located at each site.
address the many hypotheses embodied in the meta-model. Individual plots with Certain variables were constructed within the structural equation modeling
greater than 10% woody plant cover were omitted from consideration to maintain process using the composite index development methods of ref. 41. Consideration
comparability in total biomass across plots. This step resulted in the removal of of the ideas conveyed by the meta-model (Extended Data Fig. 1) and the specific
73 plots, leaving 1,126 plots in the data set analysed. Four plots were omitted owing situation being modelled suggested the need to develop index variables for soil
to incomplete plant data and one for incomplete light data. For two of the sites, live fertility and soil suitability. Soil fertility indices were developed using all meas-
mass was estimated from total mass using available information on the proportion ured soil properties and were operationally defined as the drivers of productivity,
of live to total. One apparent measurement error was detected for light data and controlling for all other effects on productivity in the model. Two indices were
the associated plot removed from the analysed sample. Random imputation meth- developed, one for site-to-site variations and another for plot-to-plot variations.
ods34 were used for cases where there were missing soil measurements at a site. Similarly, soil suitability indices were developed for the site- and plot-level data
The decision to use this approach was based on weighing the demerits of deleting using all measured soil properties as potential contributors and operationally
nearly complete multivariate data records versus introducing a modest amount of defined as the drivers of richness, controlling for all other effects on richness in
random error through the imputation process. the model.
Study plots in this investigation had a perimeter of 5 m × 5 m and were separated Modelling with composites in structural equation models involved a two-step
by 1 m walkways. A single 1 m × 1 m subplot within each plot was permanently process. First, we constructed a fully specified structural equation model (as repre-
marked and sampled for species richness during the season of peak biomass. Sites sented in Fig. 2), but providing a specific set of soil properties to serve as formative
with strong seasonal variation in composition were sampled twice during the indicators for soil fertility and soil suitability. Variables that did not contribute to
season to assemble a complete list of species. To obtain an estimate of site-level the total model (on the basis of model fit indices) were eliminated individually for
richness, we used a jack-knife procedure35. (Because there have been some recent the two composites being formed. The resulting prediction equations were used to
advances in the reduction of certain sources of bias in richness estimation36, we compute index scores. Then, the model was reconstructed, substituting the indi-
checked our original results by computing site-level richness using the new iNEXT ces in place of the collection of individual soil properties. Documentation of this
R package. The correlation between the two estimates of richness was found to process is provided in the Supplementary Information computer code (R script).
be 0.972.) Analyses. A structural equation model was developed based on the ideas embod-
Productivity and total above-ground biomass were sampled immediately adja- ied in the meta-model, available data, and the principles and procedures laid out
cent to the permanent vegetation subplot. Vegetation was sampled destructively in ref. 42. Indicators for constructs were chosen from the set of variables availa-
by clipping at ground level all above-ground biomass of individual plants rooted ble and quantities that could be computed from them (Extended Data Table 1).
within two 0.1-m2 (10 cm × 100 cm) strips. Harvested plant material was sorted The modelling approach used was semi-exploratory in that while we worked to
into the current year’s live and recently senescent material, and into previous year’s address the general hypothesis embodied in the meta-model, the precise variables
growth (including litter). For shrubs and sub-shrubs, the current year’s leaves and (for example, mean annual precipitation versus mean annual precipitation in the
stems were collected. Plant material was dried at 60 °C to a constant mass and warmest quarter of the year) to use for certain constructs (specifically, resource
weighed to the nearest 0.01 g. We used the current year’s biomass increment as supplies and regulators) were determined empirically. Compositing techniques
our estimate of annual above-ground productivity, which commonly serves as were used to estimate construct-level effects41. For comparative purposes, we ana-
a measurable surrogate for total productivity37,38. All sites used this protocol to lysed the bivariate pattern in Fig. 1A using a variety of regression models, including
estimate productivity (except for the Sevilleta, New Mexico, site which relied Ricker-type nonlinear models as well as second- and third-order polynomials.
on species-specific allometric relationships39). Total above-ground biomass was A three-parameter Ricker-type model provided the best fit for the data.
computed as the sum of the current year’s biomass and that from previous years Data were screened for distributional properties and nonlinear relations.
and included remaining dead material (litter). Photosynthetically active radiation Several variables were log-transformed as a result of evaluations (Extended Data
was measured at the time of peak biomass, both above the vegetation and at the Table 1). We used the R software platform43 and the lavaan package44 along
ground surface, the ratio representing the proportion of available light reaching with the lavaan.survey45 package for our structural equation model analyses.
the ground. Degree of shading was computed as 1.0 minus the proportion of light For the plot-scale model, robust χ2 tests, as implemented in the [Link]
reaching the ground. package, were used to judge variable inclusion and model adequacy because
Within each plot, 250 g of soil were collected and air dried for processing of the nested nature of the plot-level data. Each link in the final model was
and soil archiving. Total soil %C and %N were measured using dry combustion evaluated for significant contribution to the model. Final model fit to data was
gas chromatography analysis (COSTECH ESC 4010 Element Analyzer) at the very good for both submodels. Model fit indices were supplemented by using
University of Nebraska. All other soil analyses were performed at A&L Analytical additional diagnostic evaluations that involve visualizing residual relationships

© 2016 Macmillan Publishers Limited. All rights reserved

 RESEARCH Letter

to evaluate conditional independence29. These residual visualizations allowed, Code availability. The computer script associated with the analyses in this paper
among other things, an ability to evaluate linearity assumptions and implement is available as part of the Supplementary Information.
curve-fitting procedures if needed (which was only the case for the composite
relationships in this case). Our structural equation model in this case is non- 31. Grace, J. B., Anderson, T. M., Olff, H. & Scheiner, S. M. On the specification of
recursive and includes a causal loop. Models of this form are commonplace in structural equation models for ecological systems. Ecol. Monogr. 80, 67–87
(2010).
structural equation model applications, although they come with some additional 32. Regan, H. M., Colyvan, M. & Burgman, M. A. A taxonomy and treatment of
assumptions and requirements. Specifically, there is a requirement for unique uncertainty for ecology and conservation biology. Ecol. Appl. 12, 618–628
predictors for the elements involved in loops, a requirement that was met in this (2002).
case. Additional analysis details are documented in the R script used for the 33. Borer, E. T. et al. Finding generality in ecology: a model for globally distributed
analysis (Supplementary Information). experiments. Methods Ecol. Evol. 5, 65–73 (2014).
34. Gelman, A. & Hill, J. Data Analysis Using Regression and Multilevel/Hierarchical
Multi-level relations were incorporated into the architecture of our model.
Models (Cambridge Univ. Press, 2007).
Several ways to incorporate both site- and plot-level variations in the model were 35. Chao, A. Estimating the population size for capture-recapture data with
considered and multiple approaches evaluated to ensure results are general. In the unequal catchability. Biometrics 43, 783–791 (1987).
model form presented, we chose to follow modern hierarchical modelling prin- 36. Chao, A. & Jost, L. Coverage-based rarefaction and extrapolation:
ciples and allow plot-level observations to depend on site-level parameters, since standardizing samples by completeness rather than size. Ecology 93,
plots were nested within sites. The result of choosing this approach means site-level 2533–2547 (2012).
37. Lauenroth, W., Hunt, H., Swift, D. & Singh, J. Estimating aboveground net
explanatory effects can filter down to the plot level while plot-level explanatory primary production in grasslands: a simulation approach. Ecol. Modell. 33,
variables (for example, pathways from edaphic conditions to plot richness) explain 297–314 (1986).
additional plot-to-plot variations in responses that are not predicted from site- 38. Oesterheld, M. & McNaughton, S. J. in Methods in Ecosystem Science (eds
level (mean) conditions. Consistent with the capabilities of the structural equation Sala, O. E., Jackson, R. B., Mooney, H. A. & Howarth, R. W.) Ch. 10, 151–157
model software used in our analyses (described below), we estimated site- and (Springer, 2000).
39. Muldavin, E. H., Moore, D. I., Collins, S. L., Wetherill, K. R. & Lightfoot, D. C.
plot-level submodels using a two-stage approach, first estimating parameters for Aboveground net primary production dynamics in a northern Chihuahuan
the site-level component and then using site productivity, biomass, and richness as Desert ecosystem. Oecologia 155, 123–132 (2008).
exogenous predictors in the plot-level component. Comparisons with results from 40. Hijmans, R., Cameron, S., Parra, J., Jones, P. & Jarvis, A. WorldClim, version 1.3
separate site- and plot-level models led to very similar conclusions, although the (Univ. California, Berkeley, 2005).
hierarchical approach used allowed a better integration of processes and greater 41. Grace, J. B. & Bollen, K. A. Representing general theoretical concepts in
structural equation models: the role of composite variables. Environ. Ecol. Stat.
variance explanation.
15, 191–213 (2008).
One of our objectives in this study was to assess the model dimensionality 42. Grace, J. B., Scheiner, S. M. & Schoolmaster, D. R. Jr in Ecological Statistics (eds
needed to detect the hypothesized signals in the data. To do this, we started with Fox, G. A., Negrete-Yankelevich, S. & Sosa, V. J.) Ch. 8, 168–199 (Oxford Univ.
the most complete model (Fig. 2) and eliminated variables from the model (always Press, 2015).
retaining richness and some measure of biomass production, either productivity 43. R Development Core Team. R: a language and environment for statistical
or total biomass). We then made any modifications needed to ensure adequate computing (R Foundation for Statistical Computing, 2012).
44. Rosseel, Y., Oberski, D., Byrnes, J., Vanbrabant, L. & Savalei, V. lavaan: latent
model-data fit for these reduced-form models. The consequences of model simpli- variable analysis (software) (R Foundation for Statistical Computing, 2013).
fication was judged on the basis of signal retention, in particular a loss of capacity 45. Oberski, D. [Link]: an R package for complex survey analysis of
to detect signals associated with the remaining parts of the model. structural equation models. J. Stat. Softw. 57, 1–27 (2014).

© 2016 Macmillan Publishers Limited. All rights reserved

 Letter RESEARCH

Extended Data Figure 1 | Structural equation meta-model showing Literature and meta-model development are discussed in the
hypothesized probabilistic expectations based on literature related Supplementary Information. Specific implementations of this generalized
to the productivity–diversity debate. Solid lines represent expected model for particular cases will probably differ in detail as appropriate for
positive effects, dashed lines represent expected negative effects. the situation and available data.

© 2016 Macmillan Publishers Limited. All rights reserved

 RESEARCH Letter

Extended Data Table 1 | Model variables and their indicators*



0RGHO9DULDEOHV ,QGLFDWRU9DULDEOHV 8QLWV 

6LWH5LFKQHVVORJ HVWLPDWHGVLWHULFKQHVV QXPEHUVLWH


6LWH%LRPDVVORJ SHDNVHDVRQWRWDODERYHJURXQGELRPDVVLQFOXGLQJOLWWHU PHDQJP 

6LWH3URGXFWLYLW\ORJ SHDNVHDVRQOLYHDERYHJURXQGELRPDVVLQFUHPHQW PHDQJP 

3ORW5LFKQHVVORJ VSHFLHVLQDSORW QXPEHUSORW

6KDGLQJORJ SURSRUWLRQDOUHGXFWLRQLQOLJKWDWJURXQGVXUIDFH ±

3ORW%LRPDVVORJ SHDNVHDVRQWRWDODERYHJURXQGELRPDVVLQFOXGLQJOLWWHU JPSHUSORW

3ORW3URGXFWLYLW\ORJ SHDNVHDVRQOLYHDERYHJURXQGELRPDVVLQFUHPHQW JPSHUSORW

&OLPDWH(IIHFWRQ6LWH5LFKQHVV PHDQSUHFLSLWDWLRQLQZDUPHVWTXDUWHU PP

&OLPDWH(IIHFWRQ6LWH3URGXFWLYLW\ PHDQSUHFLSLWDWLRQLQZDUPHVWTXDUWHUWHPSHUDWXUHLQZHWWHVWTXDUWHU PPR&

'LVWXUEDQFH(IIHFWRQ6LWH5LFKQHVV KLVWRU\RIPDMRUDQWKURSRJHQLFLQIOXHQFHV RU

'LVWXUEDQFH(IIHFWRQ6LWH%LRPDVV KHUELYRU\EDVHGRQH[FORVXUHVWXGLHV SURSRUWLRQ

+HWHURJHQHLW\ &9IRUYDULDWLRQVLQYHJHWDWLRQGHQVLW\DPRQJSORWVDWDVLWH XQLWOHVV

H[SUHVVHGLQWHUPVRIFDQRS\OLJKWLQWHUFHSWLRQYDULDWLRQV

6RLO6XLWDELOLW\VLWHDQGSORWOHYHO IXQFWLRQRIVRLO13&WH[WXUHS+WKDWPD[LPL]HVULFKQHVV XQLWOHVV

FRQWUROOLQJIRURWKHUFRQGLWLRQV

6RLO)HUWLOLW\VLWHDQGSORWOHYHO IXQFWLRQRIVRLO13&WH[WXUHS+WKDWPD[LPL]HVSURGXFWLYLW\ XQLWOHVV

FRQWUROOLQJIRURWKHUFRQGLWLRQV
*The data are provided along with the Supplementary Information.

**Units given are for the raw (untransformed) variables.




© 2016 Macmillan Publishers Limited. All rights reserved

 Letter RESEARCH

Extended Data Table 2 | Results of model dimensionality evaluations

Models of different complexity were evaluated to determine the potential for model simplification. The bases for comparison were ‘full models’ for each site, as shown in Fig. 2. The consequences of
removing various components of the models are summarized under the columns ‘Signals lost’ and ‘R2 for richness’.

© 2016 Macmillan Publishers Limited. All rights reserved

 RESEARCH Letter

Extended Data Table 3 | Basic information on the study sites included in the final analyses

A total of 39 sites from the Nutrient Network ([Link] possessed sufficiently complete multivariate data to be incorporated into this analysis.

© 2016 Macmillan Publishers Limited. All rights reserved