1
Microbial functional diversity enhances predictive models linking environmental parameters to
2
ecosystem properties
3 4
Jeff R. Powell1*, Allana Welsh2, Sara Hallin3
5 6 7
1
8
Australia
9
2
Department of Geology, University of Illinois at Urbana Champaign, Urbana, IL, USA
10
3
Swedish University of Agricultural Sciences, Department of Microbiology, Uppsala, Sweden
Hawkesbury Institute for the Environment, University of Western Sydney, Richmond,
11 12
* Author for correspondence: +61 (0)2 4570 1093,
[email protected]
13 14 15
Running title: Bacterial diversity drives nitrogen cycling
16 17
Manuscript type: Article
18 19 20
Words (abstract): 146
21
Words (main text): 3798
22
References: 48
23
Figures: 4
24
Supplemental Tables: 1
25
26
Abstract
27
Microorganisms drive biogeochemical processes, but linking these processes to real changes in
28
microbial communities under field conditions is not trivial. Here, we present a model-based
29
approach to estimate independent contributions of microbial community shifts to ecosystem
30
properties. The approach was tested empirically, using denitrification potential as our model
31
process, in a spatial survey of arable land encompassing a range of edaphic conditions and two
32
agricultural production systems. Soil nitrate was the most important single predictor of
33
denitrification potential (ΔAICc = 20.29); however, the inclusion of biotic variables
34
(particularly the evenness and size of denitrifier communities [ΔAICc = 12.02], and the
35
abundance of one denitrifier genotype [ΔAICc = 18.04]) had a substantial effect on model
36
precision, comparable to the inclusion of abiotic variables (R2biotic = 0.28, R2abiotic = 0.50, R2biotic
37
+ abiotic
38
communities to ecosystem functioning. By making this link, we have demonstrated that
39
including aspects of microbial community structure and diversity in biogeochemical models can
40
improve predictions of nutrient cycling in ecosystems and enhance our understanding of
41
ecosystem functionality.
= 0.76). This approach provides a valuable tool for explicitly linking microbial
42 43
keywords: ecosystem services, functional diversity, model selection, multimodel inference,
44
nitrogen cycling
45 46 47 48
2
49
Introduction
50 51
Microorganisms are dominant components in the functioning of all ecological systems, and
52
their roles in nutrient cycling are essential to carbon dynamics and the maintenance of
53
ecosystem productivity (van der Heijden et al. 2008). As such, some ecosystem models of
54
carbon and nutrient cycling include parameters representing microbial abundance (e.g.,
55
microbial biomass, bacterial:fungal ratios, or coarse functional guilds), but these parameters
56
only capture variation at very coarse scales (Treseder et al. 2013). Parameters representing
57
microbial diversity, however, are entirely absent from predictive models of ecosystem
58
functioning, even though recent work demonstrates that inclusion of microbial processes in
59
these models can improve projections (Allison 2012, Wieder et al. 2013, Kaiser et al. 2014,
60
Reed et al. 2014). Moreover, the composition and diversity of microbial communities has had
61
important consequences on ecosystem process rates in manipulative experiments (Bell et al.
62
2005, van der Heijden et al. 1998, Maherali and Klironomos 2007, Wittebolle et al., 2009).
63 64
Experimental demonstrations of microbial community contributions to ecosystem functioning
65
have been conducted under varying levels of complexity, from laboratory microcosms to field
66
mesocosms, but usually in a context in which the initial microbial communities have been
67
eliminated or highly disturbed in order to establish assemblages of known complexity and
68
which may be dominated by specific taxa. In many, if not all cases, researchers are measuring
69
ecosystem properties during early stages of assembly, which may lead to biases in the effect
70
sizes associated with diversity treatments (Fukami et al. 2010, Dickie et al. 2012). In addition,
71
establishment of experimental treatments is performed either using a small subset of microbial
72
taxa, representing the minority that are actually cultivable from environmental samples (e.g.,
73
Bell et al. 2005), or using black box approaches, like serial dilutions of soils, that are highly
3
74
effective at reducing diversity but potentially confound diversity with other factors such as
75
species identity and evenness (e.g., Garland and Lehman 1999, Philippot et al. 2013, Wertz et
76
al. 2006). To confirm the relevance of these results and to demonstrate that the interdependence
77
of microbial community diversity and functioning are not artifactual, we require studies of
78
communities that have assembled under those environmental conditions in which microbial
79
functional diversity is hypothesised to be important (Hallin et al. 2012, Maherali and
80
Klironomos 2007, Pfisterer et al. 2004).
81 82
The most common approach to study the functional consequences of microbial diversity in the
83
environment is to correlate structural responses in microbial communities to environmental
84
parameters and functional properties of divergent ecosystem types, using bivariate correlations
85
or multivariate ordinations (Torsvik and Øvreås 2002). There is a certain degree of circularity
86
in this approach as aspects of the ecosystem drive shifts in microbial communities, and these
87
structured communities are then suggested to have functional consequences within the
88
ecosystem (such as for plant-soil feedbacks [Bever et al. 1997] or for biogeochemical process
89
rates [Beare et al. 1995]). Thus, the causal nature of these interactions is drawn into question,
90
particularly when variation in microbial communities is imposed via experimental
91
manipulations or other anthropogenic forces that are likely to influence other biotic and abiotic
92
components of the ecosystem. Plant ecologists have faced the same issues while attempting to
93
incorporate plant functional diversity attributes into ecosystem service assessments (e.g. Loreau
94
et al. 2001). Functional diversity of plant communities has been demonstrated to be an integral
95
component of ecosystem functioning in controlled experiments (beginning with Tilman et al.
96
1997 and in many subsequent experiments), but quantifying these effects under field conditions
97
in natural systems is a challenge (Romanuk et al. 2009). One approach, developed by Diaz et al.
98
(2007), utilises multiple regression and involves two stages: i) identifying individual variables
4
99
associated with the abiotic environment and functional diversity that are statistically significant
100
predictors of an environmental property in generalized linear models, and ii) combining these
101
variables into a single predictive model and inferring the contributions of each variable from
102
their estimated coefficients in these models. While still fundamentally correlative, this approach
103
gets closer to demonstrating 'causality' as used by Shipley (2002), i.e. fitting data to models
104
representing causal hypotheses and comparing this fit to that obtained using models
105
representing appropriate alternative, mechanistic hypotheses.
106 107
Here, we extend the approach developed by Diaz et al. (2007) to attributes of microbial
108
functional diversity, including community-level mean trait values (e.g., functional gene copy
109
number representing genetic potential of the specific function), trait distributions (implicitly
110
linked to estimates of community richness and evenness), and idiosyncratic species effects
111
(based on patterns of genotype presence/absence or relative abundance). We also improve on
112
the statistical approach by incorporating multimodel-inference into the procedure, which
113
explicitly estimates variable importance and determines the sensitivity of coefficient estimates
114
to model structure (Burnham and Anderson, 2002). We then include an empirical test of this
115
approach, estimating the contribution of functional diversity in denitrifying communities to
116
variation in denitrification potential across 44 ha of arable land encompassing a range of
117
edaphic conditions and two agricultural production systems (Enwall et al. 2010). The original
118
publication observed several correlations among denitrification potential and the abundance and
119
composition of denitrifier communities, but also observed significant covariation between these
120
properties and several descriptors of the abiotic environment, representing a clear example of
121
the challenges of interpreting correlations from multivariate data collected in an observational
122
context.
123
5
124 125
Materials and Methods
126 127
Model structure
128 129
We started with the multiple regression procedure developed by Diaz et al. (2007), a six-step
130
procedure based on the mass ratio hypothesis linking ecosystem functioning to the traits of
131
dominant (in terms of contributions to biomass) plant species (Grime 1998). Briefly, we
132
classified predictor variables into one of four classes: i) abiotic environmental variables, ii)
133
variables representing aggregate trait values for the community (representing the community-
134
weighted mean trait value), iii) variables representing the distribution of trait values for the
135
community (representing functional diversity), and iv) variables representing the relative
136
abundance or presence-absence of specific genotypes in the community (representing
137
idiosyncratic species effects). Justifications for the classification of the specific variables used
138
here are provided below (under 'Empirical example'). For the first three classes, we included all
139
relevant predictor variables in a linear model explaining variation in the response variable. For
140
the fourth class, we employed pairwise correlations between the response variable and each
141
predictor variable (Gotelli et al. 2011) employing the Šidák correction for multiple
142
comparisions (Šidák 1967) because the number of predictors (genotypes) greatly outnumbered
143
the number of samples. For Diaz’ fifth step, those predictor variables that were significant in
144
each of the previous four steps were included in another linear model, accounting for the
145
combined effects of abiotic factors and functional diversity, that explained variation in the
146
response variable and the significance of each predictor assessed. The sixth step proposed by
147
Diaz et al. (2007) accounts for nonlinear relationships between the response and predictor
6
148
variables. We included this step, but did not find support for nonlinearities in these data, so we
149
do not discuss this step further.
150 151
To test the significance of predictor variables, Diaz et al. (2007) used an ascending procedure,
152
selecting the most parsimonious model in the final step based on the Akaike Information
153
Criterion (AIC). Instead of basing our inferences on this single model, we used multimodel
154
inference and model averaging to assess the contribution of each variable, and class of
155
variables, to efficiently explain variation in the response variable. In our procedure, predictors
156
were selected randomly for inclusion in each model, model fit was estimated using AIC,
157
corrected for small sample size (AICc; Hurvich & Tsai 1989), then model-averaged estimates
158
of standardised effect sizes and variable 'importance' were obtained across all models within
159
each subset using Akaike weights (Burnham and Anderson 2002). Standardised effect sizes of
160
predictor variables were estimated from the weighted averages of t-statistics. The 'importance'
161
of each predictor variable was estimated by calculating the sum of Akaike weights for models
162
in which they were included; individual weights are expressed as a proportion, with the weights
163
across all models in the set summing to one. We used an 'importance' threshold of 0.5 to justify
164
inclusion of predictor variables in subsequent models. We preferred this procedure over the
165
ascending procedure used by Diaz et al. because it allowed us to account for situations in which
166
i) the observation of effects associated with one variable was conditional on the inclusion of
167
another variable in the model and ii) there were not enough degrees of freedom in the model to
168
account for all predictor variables. We ranked the importance of each class of variables based
169
on the estimated effect on Akaike weights when all relevant variables were dropped from the
170
full model.
171 172
Empirical example
7
173 174
To demonstrate the utility of the approach described above, we used a dataset derived from a
175
spatial survey of two agricultural production systems (integrated management: optimising
176
nutrient inputs and pesticide use; organic management: no pesticide or mineral fertiliser use,
177
with nitrogen inputs by cultivation of legumes) within the Logården research farm in Sweden
178
(58°20'N, 12°38'E; altitude, 50 m) and previously used to infer soil factors driving spatial
179
patterns in denitrifier communities (Enwall et al. 2010). Each farming system consists of a
180
seven-year crop rotation spread across seven fields. We focussed on potential denitrification
181
rate (estimated using the acetylene inhibition technique of Enwall et al. [2005]) as the
182
ecosystem response. The rates do not reflect the actual rates, but rather the enzymatic potential
183
in the system since the assay is conducted under optimal conditions, with supplemental
184
substrate in excess and under anaerobic conditions. Data were collected at 56 sampling
185
locations within the farm and specific details on sample and data collection can be found in the
186
original publication (Enwall et al. 2010). Briefly, measurements included a variety of soil
187
chemical (total-N, nitrate-N, ammonia-N, and dissolved organic N; total organic C and
188
dissolved organic C; plant-available P, K, Mg, and Ca; total Cu; pH) and physical (soil
189
penetration resistance at 6, 8, 14, and 16 cm depth; water holding capacity at 0.5 and 5 kPa;
190
water infiltration rate; pore volume; clay content; bulk density; water content at sampling)
191
characteristics. Microbiological parameters included the genetic potential of denitrification in
192
terms of size and structure of the denitrifier communities based on DNA extracted from the
193
soil. These parameters were determined using qPCR and terminal restriction fragment length
194
polymorphism (T-RFLP), respectively, of the genes nirS and nirK encoding the two types of
195
nitrite reductases frequently used as functional markers for denitrifiers (Throbäck et al, 2004).
196
Denitrification is an anaerobic respiration pathway and the defining step is the reduction of
197
nitrite to nitric oxide catalyzed by one of two different types of nitrite reductase, NirK or NirS,
8
198
which are mutually exclusive in denitrifying organisms with a few exceptions (Graf et al.
199
2014). When necessary, variables were log10-transformed to ensure normality of error residuals.
200 201
For each of the nirS and nirK communities, richness and Pielou's evenness (Shannon diversity /
202
log richness; Pielou 1966) of nirS and nirK genotypes, based on the number and relative peak
203
areas of terminal restriction fragments from the TRFLP-analysis, were used as a proxy for
204
functional trait diversity. This is a reasonable assumption, given that functional diversity in
205
bacteria has been observed at a range of taxonomic levels that are likely to be captured by
206
fingerprinting approaches (e.g., Torsvik and Øvreås 2002, Mendes et al. 2011). In the absence
207
of actual trait data from members of each genotype observed in environmental surveys, an
208
unrealistic (if not impossible) expectation using available technology, using structural
209
descriptions of microbial communities to represent functional diversity is a reasonable way
210
forward. Idiosyncratic species effects were modelled by regressing potential denitrification
211
rates against the relative abundance of genotypes; we used the Šidák correction for multiple
212
comparisions (Šidák 1967) because of the large number of genotypes (nirK: 113; nirS: 68). In
213
cases where a specific genotype was not observed in multiple samples, we also performed a t-
214
test to compare denitrification potential in samples where a genotype was present to those
215
where it was absent. This latter step was performed with the intention to exclude genotypes
216
where a significant relationship with relative abundance was observed but model fit was poor;
217
however, we did not encounter this scenario so the step was not necessary for these data. Due to
218
differences in the number of variables present in each subset, model-averaged parameter
219
estimates were generated from a variable number of models at each step of the procedure. In the
220
first step, up to twelve abiotic variables were randomly selected for each model and parameters
221
estimated across 16102 models. In the second through the fifth step, all potential combinations
222
of parameters were tested resulting in four, 64, 16, and 16384 models, respectively. Once the
9
223
best predictive model was identified, we also compared the fit of this model to the data with an
224
equivalent model that also included a term representing the crop production system (integrated
225
management or organic management) associated with the field from which each sample was
226
collected. In addition, we estimated the decrease in model fit after dropping individual terms
227
from the model. Finally, we examined bivariate correlations between the response variable and
228
each predictor variable in the best predictive model. All analyses were performed using R (R
229
Core Development Team 2013); the ‘vegan’ package (Oksanen et al. 2013) was used to
230
estimate diversity metrics.
231 232 233
Results
234 235
We found evidence for all four classes of variables (abiotic environmental variables, aggregated
236
trait values within communities, functional diversity, and idiosyncratic species effects) being
237
important predictors of potential denitrification rates. Initial independent screening of the four
238
classes of variable identified several variables that were not deemed important when the
239
combined effects across the four variables were assessed (Figure 1). This suggests that these
240
variables may have had indirect effects on denitrification rates via other variables in the model
241
or that they may have been correlated to denitrification due to an independent response to
242
variation in another variable.
243 244
Of the important predictors, nitrate concentration was the most important single predictor of
245
denitrification potential (ΔAICc = 20.29; Figure 2). The inclusion of abiotic variables had a
246
greater effect on the explanatory power of the model than the inclusion of biotic variables
247
(ΔAICc = 44.91 vs. ΔAICc = 26.57, respectively; Figure 2). However, the effects of including
10
248
the biotic variables as predictors was still substantial, accounting for approximately one-third of
249
the variation explained by the full model (R2 = 0.76; Figure 2). Of these, the evenness of nirS
250
communities and the abundance of one nirS genotype had a comparable level of importance
251
(ΔAICc = 12.02 and 18.04, respectively) to the abiotic variables, while the abundance of nirK
252
gene copies was less important (ΔAICc = 3.19). The bivariate correlations between most of
253
these variables and denitrification rate suggested positive, linear relationships (Figure 3). The
254
exceptions to this were the predictors relating to soil penetration resistance; these had a more
255
complex relationship with denitrification, with vertical structure explaining variation. Rates
256
showed a very weak, non-significant tendency to decrease with increasing penetration
257
resistance at 6 cm depth but increase with penetration resistance at 8 cm depth (Figure 3).
258
However, we observed a positive relationship between denitrification and the difference in
259
penetration resistance at these two depths (R2 = 0.15; Figure 4).
260 261
Adding a term that accounted for the crop production system (integrated versus organic
262
management) did not improve model fit and, in fact, resulted in a large reduction in model fit
263
relative to the best model in Figure 2 according to the information criterion (ΔAICc = 2.82).
264
Inspection of the model coefficients confirmed that each of the biotic and abiotic variables
265
remained statistically significant predictors of variation in potential denitrification rate, but that
266
the crop production system was not a significant predictor (Table A1).
267 268 269
Discussion
270 271
Parameters representing microbial diversity are seldom considered in predictive models of
272
ecosystem functioning. With an improved model-based approach and using denitrification
11
273
potential and the microbial community corresponding to this function as a model system, we
274
provide a direct demonstration that variation in microbial community structure under field
275
conditions can have functional consequences at ecologically relevant spatial scales; a result that
276
has only been shown previously using indirect approaches and under highly controlled
277
conditions (e.g., van der Heijden et al. 1998, Wittebolle et al., 2009). Environmental parameters
278
and microbial community structure (evenness of nirS communities) both predicted significant,
279
substantial, and nearly similar amounts of variation in the functional response. Our results show
280
that ignoring variation in microbial community structure limits our ability to both predict and
281
understand ecosystem functionality as the biotic parameters accounted for about one-third of
282
the variation explained by the model. The approach also served to separate variables that were
283
predictive of denitrification potential from others that were collinear. While significant
284
correlations among these variables were observed in the original analysis of these data (Enwall
285
et al. 2010), many more correlations were also observed to be significant that our current
286
analysis identified as not being directly associated with potential denitrification rates.
287 288
The nirS-denitrifier community evenness and a specific nirS genotype were the most important
289
biotic variables predicting function, while the importance of the size of the nirK-denitrifier
290
community was somewhat reduced in comparison. This agrees with what was observed by
291
correlations in the original data set (Enwall et al. 2010) and in other studies reporting either a
292
significant and/or stronger correlation between denitrification activity and nirS abundance when
293
compared to nirK (Philippot et al. 2009; Graham et al. 2010). It has been suggested that nirS-
294
and nirK-type denitrifier communities occupy different ecological niches (e.g., Santoro et al.,
295
2006; Desnues et al., 2007; Smith & Ogram, 2008; Yuan et al., 2012) and coexisting nirS and
296
nirK communities are likely not under the same community assembly rules (Jones and Hallin,
297
2010). In another experimental agricultural system, Hallin et al. (2012) observed that
12
298
denitrification process rates could be maintained over a range of temperatures and salinity
299
levels and this functional operating range was expanded for bacterial communities that
300
assembled under certain soil management practices and contained key genotypes. Others have
301
also observed key bacterial genotypes to be correlated with process rates in communities
302
assembling under field conditions (such as methanotrophy, Nazaries et al. 2011), but these
303
observational approaches do not account for co-linearity with environmental parameters that
304
may also drive these processes (e.g., nitrate concentration is a primary and direct determinant of
305
soil processes such as denitrification [observed here] and methane cycling [Reay et al. 2001]).
306
For instance, four nirS genotypes were observed to be significantly correlated with potential
307
denitrification rate, but these all covaried in abundance (Appendix B) and only two were
308
determined to explain significant variation when assessed together. The impacts on functioning
309
of key denitrifying genotypes may be related to enhanced resource extraction efficiencies or
310
successful interactions with that genotype (Salles et al. 2009). By using our model based
311
approach, we could demonstrate that such changes in the abundance of key genotypes can have
312
direct effects on soil processes, independent from changes in edaphic properties.
313 314
We observed that the evenness of denitrifier communities (particularly nirS communities) was
315
more strongly associated with potential denitrification rates than was denitrifier richness. Some
316
researchers have found denitrifier richness to have low explanatory power when predicting
317
denitrification; functional dissimilarity among genotypes was a better predictor than richness in
318
experimental assemblages of up to 16 genotypes (Salles et al. 2009). Others have observed
319
positive effects associated with denitrifier diversity, although the extent that this is due to
320
changes in richness and evenness is not clear. Philippot et al. (2013) used a black-box approach
321
(serial dilutions of intact microbial communities) to achieve a four- to five-fold reduction in
322
denitrification rate, along with a 75 % decrease in denitrifier richness and a 30 % decrease in
13
323
evenness. Based on visual comparisons of the estimates across three levels of dilution, the
324
roughly linear decrease in denitrification rate appeared to track better with the roughly linear
325
decrease in richness than the nonlinear decrease in evenness, but this is speculative as the
326
relationship was not tested in an explicit fashion. Similar to what we report here, community
327
evenness has been previously linked to positive effects on ecosystem functioning, but the
328
mechanism by which it can confer these effects is poorly understood. In one example, bacterial
329
communities exhibited greater functionality in terms of net ecosystem denitrification with
330
increasing evenness, particularly under selective stress (Wittebolle et al. 2009). This result may
331
have been due to bacterial communities with high evenness being more likely to contain
332
sufficient stress-tolerant populations capable of responding to stress conditions, analogous to
333
selection effects that arise due to dominance by species with particular traits that directly impact
334
ecosystem processes (Loreau and Hector 2001). Selection effects can also be dependent on a
335
particular set of environmental conditions, which may be ephemeral; thus, communities that are
336
co-dominated by multiple species, each contributing to ecosystem functioning but under
337
different circumstances, are better able to maintain function (Isbell et al. 2011). Relationships
338
between evenness and ecosystem functions may also emerge from interactions among
339
community members when dominance is reduced. Dimitrakopoulos (2010) observed that low
340
litter-species evenness in mixtures containing litter from multiple herbaceous plant species was
341
linked to reductions in the rate of litter decomposition, partially explained by the contributions
342
of plant community composition to C:N ratios. Our observations do not tell us why we
343
observed a positive relationship between denitrification rates and evenness, independently of
344
community composition. It is possible that a selective environmental stress was present at the
345
site during or prior to sampling but not accounted for in the measured parameters, and that
346
variation in bacterial community structure and composition mediated the functional response
347
(Hallin et al. 2012).
14
348 349
The importance of appropriate model selection approaches and multimodel inference was
350
highlighted when considering abiotic predictors of denitrification rates in the model,
351
particularly penetration resistance, an indicator of increasing soil compaction, at different soil
352
depths. One could expect that soil compaction would have positive relationship with
353
denitrification, being an anaerobic process. In our research, there was a complex relationship
354
between soil penetration resistance at two depths suggesting that variation in the vertical
355
structure of soils may have important consequences for nutrient cycling. While relationships
356
between potential denitrification rates and penetration resistance at 6 or 8 cm depth were weak
357
(Figure 3), denitrification rates were observed to increase substantially as penetration resistance
358
increased moving between these two layers, perhaps reflecting small-scale variation in the
359
movement and retention of water and/or nitrate in soil. A forward-stepwise approach that
360
accounted for each of these variables separately would have missed this complex relationship,
361
which explained up to 12 % of the variation in denitrification rate when included in the model
362
with the other predictors (compare R2 values across the models in Figure 2). In addition,
363
inferring effect size estimates associated with penetration resistance predictors across multiple
364
models also revealed large contributions of these variables to variation in denitrification rates
365
(Figure 1). Such a multimodel approach allows the user to observe whether the direction and
366
magnitude of the relationships between response variables and predictors changes depending on
367
the inclusion of other variables in the model, resulting in estimates that are closer to zero than
368
in the single best model (Burnham and Anderson 2002).
369 370
Attempts to intervene in microbial ecological processes under field conditions may be informed
371
by demonstrating the relative contributions of microbial functional diversity and environmental
372
parameters to these processes. For example, the specific case presented here suggests that
15
373
altering vertical gradients in the physical properties of the soil (for instance, using tillage
374
strategies [Linn and Doran 1984]) might be used to establish the conditions in which microbial
375
communities contribute to these processes. In addition, the idiosyncratic species component of
376
the predictive models helps to identify taxa that may be targeted for management. However, if
377
the contribution of microbial community evenness to ecosystem functioning, as observed here,
378
is a general phenomenon of many microbial ecological processes, this would suggest that the
379
relative abundances of dominant and subordinate members need to be accounted for when
380
managing bacterial communities and simple inoculation strategies may not achieve any gains.
381
This latter point is a challenging proposition and requires a much greater understanding of the
382
drivers of evenness in bacterial communities.
383 384
While using this model-based approach highlights the linkages between bacterial community
385
structure and ecosystem functioning, the data used to estimate these linkages can suffer from
386
certain limitations. For instance, functional measurements often reflect ecosystem potential
387
rather than functioning. Moreover, PCR-based approaches are often criticized for being biased
388
(for example, Kanagawa 2003), and issues of resolution exist for community data generated via
389
fingerprinting methods (for example, Dickie and FitzJohn 2007). Nevertheless, the patterns
390
observed using an approach as suggested here can lead to the generation of hypotheses that,
391
when tested, will further our understanding of how microbial diversity contributes to ecosystem
392
functioning, similar to the development of knowledge surrounding the contributions of plant
393
functional diversity to ecosystem services (Diaz et al. 2007, Lavorel et al. 2011, Powell et al.
394
2013).
395 396
In conclusion, using this model-based approach designed to link functional diversity to
397
variation in ecosystem properties and to compare the importance of this diversity to other
16
398
environmental drivers, we demonstrated that variation in bacterial communities assembled
399
under field conditions had functional consequences for biogeochemical cycling. This is a novel
400
contribution that explicitly links microbial diversity to processes that are important at much
401
larger ecological scales. It also sets the stage for further quantitative insights from observational
402
data into the explicit contributions of microorganisms to ecosystem properties and services. The
403
approach can be applied to estimate the functional consequences of temporal or spatial variation
404
in microbial communities; it can also be extended to establish microbial community
405
contributions to ecosystem multifunctionality. Examples of the latter are bacterial groups
406
catalysing multiple pathways in the nitrogen cycle, such as i) bacteria capable of denitrification,
407
nitrogen fixation and dissimilatory nitrate reduction or ii) nitrifiers that are also denitrifying.
408
Eventually, these insights may lead to better predictive models of landscape-scale variation in
409
ecosystem services, by explicitly accounting for variation in microbial communities.
410 411
Acknowledgements
412 413
We thank Brajesh Singh and Christopher Jones for discussion during early stages of this work
414
and Steven Allison and two reviewers for helpful comments. This work was partly supported by
415
the Swedish Research Council Formas (contract 2013-656).
416 417 418
References
419 420 421
Allison, S. D. 2012. A trait-based approach for modelling microbial litter decomposition. Ecology Letters 15:1058–1070.
17
422
Beare, M. H., D. C. Coleman, D. A. C. Jr, P. F. Hendrix, and E. P. Odum. 1995. A hierarchical
423
approach to evaluating the significance of soil biodiversity to biogeochemical cycling.
424
Pages 5–22 in H. P. Collins, G. P. Robertson, and M. J. Klug, editors. The Significance
425
and Regulation of Soil Biodiversity. Springer Netherlands.
426
Bell, T., J. A. Newman, B. W. Silverman, S. L. Turner, and A. K. Lilley. 2005. The
427
contribution of species richness and composition to bacterial services. Nature 436:1157–
428
1160.
429
Bever, J. D., K. M. Westover, and J. Antonovics. 1997. Incorporating the soil community into
430
plant population dynamics: the utility of the feedback approach. Journal of Ecology
431
85:561–573.
432 433
Burnham, K. P., and D. R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Springer.
434
Desnues, C., V. D. Michotey, A. Wieland, C. Zhizang, A. Fourçans, R. Duran, and P. C. Bonin.
435
2007. Seasonal and diel distributions of denitrifying and bacterial communities in a
436
hypersaline microbial mat (Camargue, France). Water Research 41:3407–3419.
437
Díaz, S., S. Lavorel, F. de Bello, F. Quétier, K. Grigulis, and T. M. Robson. 2007.
438
Incorporating plant functional diversity effects in ecosystem service assessments.
439
Proceedings of the National Academy of Sciences 104:20684–20689.
440
Dickie, I. A., and R. G. FitzJohn. 2007. Using terminal restriction fragment length
441
polymorphism (T-RFLP) to identify mycorrhizal fungi: a methods review. Mycorrhiza
442
17:259–270.
443
Dickie, I. A., T. Fukami, J. P. Wilkie, R. B. Allen, and P. K. Buchanan. 2012. Do assembly
444
history effects attenuate from species to ecosystem properties? A field test with wood-
445
inhabiting fungi. Ecology Letters 15:133–141.
18
446
Dimitrakopoulos, P. G. 2010. Influence of evenness on the litter-species richness–
447
decomposition relationship in Mediterranean grasslands. Journal of Plant Ecology 3:71–
448
78.
449
Enwall, K., L. Philippot, and S. Hallin. 2005. Activity and composition of the denitrifying
450
bacterial community respond differently to long-term fertilization. Applied and
451
Environmental Microbiology 71:8335–8343.
452
Enwall, K., I. N. Throbäck, M. Stenberg, M. Söderström, and S. Hallin. 2010. Soil resources
453
influence spatial patterns of denitrifying communities at scales compatible with land
454
management. Applied and Environmental Microbiology 76:2243–2250.
455
Fukami, T., I. A. Dickie, J. Paula Wilkie, B. C. Paulus, D. Park, A. Roberts, P. K. Buchanan,
456
and R. B. Allen. 2010. Assembly history dictates ecosystem functioning: evidence from
457
wood decomposer communities. Ecology Letters 13:675–684.
458
Garland, J. L., and R. M. Lehman. 1999. Dilution/extinction of community phenotypic
459
characters to estimate relative structural diversity in mixed communities. FEMS
460
Microbiology Ecology 30:333–343.
461
Gotelli, N. J., W. Ulrich, and F. T. Maestre. 2011. Randomization tests for quantifying species
462
importance to ecosystem function. Methods in Ecology and Evolution 2:634–642.
463
Graf, D., C.M. Jones and S. Hallin, S. 2014. Intergenomic comparisons highlight modularity of
464
the denitrification pathway and underpin the importance of community structure for N2O
465
emissions. PLOS ONE in press, doi: 10.1371/journal.pone.0114118.
466
Graham, D. W., C. Trippett, W. K. Dodds, J. M. O’Brien, E. B. K. Banner, I. M. Head, M. S.
467
Smith, R. K. Yang, and C. W. Knapp. 2010. Correlations between in situ denitrification
468
activity and nir-gene abundances in pristine and impacted prairie streams. Environmental
469
Pollution 158:3225–3229.
19
470 471 472
Grime, J. P. 1998. Benefits of plant diversity to ecosystems: immediate, filter and founder effects. Journal of Ecology 86:902–910. Hallin, S., A. Welsh, J. Stenström, S. Hallet, K. Enwall, D. Bru, and L. Philippot. 2012. Soil
473
functional operating range linked to microbial biodiversity and community composition
474
using denitrifiers as model guild. PLoS ONE 7:e51962.
475
Van der Heijden, M. G. A., R. D. Bardgett, and N. M. Van Straalen. 2008. The unseen
476
majority: soil microbes as drivers of plant diversity and productivity in terrestrial
477
ecosystems. Ecology Letters 11:296–310.
478
Van der Heijden, M. G. A., J. N. Klironomos, M. Ursic, P. Moutoglis, R. Streitwolf-Engel, T.
479
Boller, A. Wiemken, and I. R. Sanders. 1998. Mycorrhizal fungal diversity determines
480
plant biodiversity, ecosystem variability and productivity. Nature 396:69–72.
481
Hurvich, C. M., and C.-L. Tsai. 1989. Regression and time series model selection in small
482 483
samples. Biometrika 76:297–307. Isbell, F., V. Calcagno, A. Hector, J. Connolly, W. S. Harpole, P. B. Reich, M. Scherer-
484
Lorenzen, B. Schmid, D. Tilman, J. van Ruijven, A. Weigelt, B. J. Wilsey, E. S. Zavaleta,
485
and M. Loreau. 2011. High plant diversity is needed to maintain ecosystem services.
486
Nature 477:199–202.
487 488 489 490 491 492
Jones, C. M., and S. Hallin. 2010. Ecological and evolutionary factors underlying global and local assembly of denitrifier communities. The ISME Journal 4:633–641. Kaiser, C., O. Franklin, U. Dieckmann, and A. Richter. 2014. Microbial community dynamics alleviate stoichiometric constraints during litter decay. Ecology Letters 17:680–690. Kanagawa, T. 2003. Bias and artifacts in multitemplate polymerase chain reactions (PCR). Journal of Bioscience and Bioengineering 96:317–323.
20
493
Lavorel, S., K. Grigulis, P. Lamarque, M.-P. Colace, D. Garden, J. Girel, G. Pellet, and R.
494
Douzet. 2011. Using plant functional traits to understand the landscape distribution of
495
multiple ecosystem services. Journal of Ecology 99:135–147.
496
Linn, D. M., and J. W. Doran. 1984. Effect of water-filled pore space on carbon dioxide and
497
nitrous oxide production in tilled and nontilled soils. Soil Science Society of America
498
Journal 48:1267.
499 500 501
Loreau, M., and A. Hector. 2001. Partitioning selection and complementarity in biodiversity experiments. Nature 412:72–76. Loreau, M., S. Naeem, P. Inchausti, J. Bengtsson, J. P. Grime, A. Hector, D. U. Hooper, M. A.
502
Huston, D. Raffaelli, B. Schmid, D. Tilman, and D. A. Wardle. 2001. Biodiversity and
503
ecosystem functioning: current knowledge and future challenges. Science 294:804–808.
504 505 506
Maherali, H., and J. N. Klironomos. 2007. Influence of phylogeny on fungal community assembly and ecosystem functioning. Science 316:1746–1748. Mendes, R., M. Kruijt, I. de Bruijn, E. Dekkers, M. van der Voort, J. H. M. Schneider, Y. M.
507
Piceno, T. Z. DeSantis, G. L. Andersen, P. A. H. M. Bakker, and J. M. Raaijmakers.
508
2011. Deciphering the rhizosphere microbiome for disease-suppressive bacteria. Science
509
332:1097–1100.
510
Nazaries, L., K. R. Tate, D. J. Ross, J. Singh, J. Dando, S. Saggar, E. M. Baggs, P. Millard, J.
511
C. Murrell, and B. K. Singh. 2011. Response of methanotrophic communities to
512
afforestation and reforestation in New Zealand. The ISME Journal 5:1832–1836.
513
Oksanen, J., F. G. Blanchet, R. Kindt, P. Legendre, P. R. Minchin, B. O’Hara, G. L. Simpson,
514
P. Solymos, M. H. H. Stevens, and H. Wagner. 2013. The vegan package. Version 2.0-10.
515
http://cran.r-project.org/web/packages/vegan/index.html
21
516
Pfisterer, A. B., J. Joshi, B. Schmid, and M. Fischer. 2004. Rapid decay of diversity-
517
productivity relationships after invasion of experimental plant communities. Basic and
518
Applied Ecology 5:5–14.
519
Philippot, L., J. Čuhel, N. P. A. Saby, D. Chèneby, A. Chroňáková, D. Bru, D. Arrouays, F.
520
Martin-Laurent, and M. Šimek. 2009. Mapping field-scale spatial patterns of size and
521
activity of the denitrifier community. Environmental Microbiology 11:1518–1526.
522
Philippot, L., A. Spor, C. Hénault, D. Bru, F. Bizouard, C. M. Jones, A. Sarr, and P.-A. Maron.
523
2013. Loss in microbial diversity affects nitrogen cycling in soil. The ISME Journal
524
7:1609–1619.
525 526 527 528 529 530 531
Pielou, E. C. 1966. The measurement of diversity in different types of biological collections. Journal of Theoretical Biology 13:131–144. Powell, J. R., I. C. Anderson, and M. C. Rillig. 2013. A new tool of the trade: plant-trait based approaches in microbial ecology. Plant and Soil 365:35–40. R Core Development Team. 2013. R: A Language and Environment for Statistical Computing. http://cran.r-project.org/ Reay, D. S., S. Radajewski, J. C. Murrell, N. McNamara, and D. B. Nedwell. 2001. Effects of
532
land-use on the activity and diversity of methane oxidizing bacteria in forest soils. Soil
533
Biology and Biochemistry 33:1613–1623.
534
Reed, D. C., C. K. Algar, J. A. Huber, and G. J. Dick. 2014. Gene-centric approach to
535
integrating environmental genomics and biogeochemical models. Proceedings of the
536
National Academy of Sciences 111:1879–1884.
537 538
Romanuk, T. N., R. J. Vogt, and J. Kolasa. 2009. Ecological realism and mechanisms by which diversity begets stability. Oikos 118:819–828.
539
Salles, J. F., F. Poly, B. Schmid, and X. L. Roux. 2009. Community niche predicts the
540
functioning of denitrifying bacterial assemblages. Ecology 90:3324–3332.
22
541
Santoro, A. E., A. B. Boehm, and C. A. Francis. 2006. Denitrifier community composition
542
along a nitrate and salinity gradient in a coastal aquifer. Applied and Environmental
543
Microbiology 72:2102–2109.
544
Shipley, B. 2002. Cause and Correlation in Biology: A User’s Guide to Path Analysis,
545
Structural Equations and Causal Inference. Cambridge University Press.
546
Šidák, Z. 1967. Rectangular confidence regions for the means of multivariate normal
547
distributions. Journal of the American Statistical Association 62:626–633.
548
Smith, J. M., and A. Ogram. 2008. Genetic and functional variation in denitrifier populations
549
along a short-term restoration chronosequence. Applied and Environmental Microbiology
550
74:5615–5620.
551
Throbäck, I. N., K. Enwall, Å. Jarvis, and S. Hallin. 2004. Reassessing PCR primers targeting
552
nirS, nirK and nosZ genes for community surveys of denitrifying bacteria with DGGE.
553
FEMS Microbiology Ecology 49:401–417.
554 555 556 557 558
Tilman, D., J. Knops, D. Wedin, P. Reich, M. Ritchie, and E. Siemann. 1997. The influence of functional diversity and composition on ecosystem processes. Science 277:1300–1302. Torsvik, V., and L. Øvreås. 2002. Microbial diversity and function in soil: from genes to ecosystems. Current Opinion in Microbiology 5:240–245. Treseder, K. K., T. C. Balser, M. A. Bradford, E. L. Brodie, E. A. Dubinsky, V. T. Eviner, K. S.
559
Hofmockel, J. T. Lennon, U. Y. Levine, B. J. MacGregor, J. Pett-Ridge, and M. P.
560
Waldrop. 2012. Integrating microbial ecology into ecosystem models: challenges and
561
priorities. Biogeochemistry 109:7–18.
562
Wertz, S., V. Degrange, J. I. Prosser, F. Poly, C. Commeaux, T. Freitag, N. Guillaumaud, and
563
X. L. Roux. 2006. Maintenance of soil functioning following erosion of microbial
564
diversity. Environmental Microbiology 8:2162–2169.
23
565 566 567
Wieder, W. R., G. B. Bonan, and S. D. Allison. 2013. Global soil carbon projections are improved by modelling microbial processes. Nature Climate Change 3:909–912. Wittebolle, L., M. Marzorati, L. Clement, A. Balloi, D. Daffonchio, K. Heylen, P. De Vos, W.
568
Verstraete, and N. Boon. 2009. Initial community evenness favours functionality under
569
selective stress. Nature 458:623–626.
570
Yuan, Q., P. Liu, and Y. Lu. 2012. Differential responses of nirK- and nirS-carrying bacteria to
571
denitrifying conditions in the anoxic rice field soil. Environmental Microbiology Reports
572
4:113–122.
573 574 575 576
Ecological Archives
577 578
Appendix A. Model output demonstrating that crop production system was not a significant
579
predictor of denitrification potential in the Logården samples and including this term did not
580
affect other terms in the model.
581 582
Appendix B. Rank-abundance distributions within denitrifier communities characterized from
583
the Logården samples and covariation in the relative abundances of abundant denitrifier
584
genotypes.
24
585
Figure Captions
586 587
Figure 1. Estimates of variable importance based on Akaike weights of models containing each
588
variable. Effect sizes are estimated from weighted averages of regression coefficients across
589
mulitple models. Variables are coloured based on belonging to one of four classes of variable
590
(red: abiotic environmental variables; blue: aggregated trait values within communities; yellow:
591
functional diversity; green: idiosyncratic species effects).
592 593
Figure 2. Summaries of model fit for the full model (that which best explains variation in
594
potential denitrification rate), models in which one term was dropped, and models in which all
595
biotic or abiotic variables were dropped. Each row represents a model containing those
596
variables for which the cells are filled; colours refer to variable classes as for Fig. 1. Models are
597
ranked by declining explanatory power (R2 and Akaike weights [wi]).
598 599
Figure 3. Bivariate correlations between variables in the full model and potential denitrification
600
rate. Colours refer to variable classes as for Fig. 1.
601 602
Figure 4. Interdependence between penetration resistance at two soil depths in the response of
603
potential denitrification rate. Potential denitrification rate was higher in samples where
604
penetration resistance (PEN) at 8 cm was higher than expected (based on the relationship
605
between penetration resistance at 6 and 8 cm).
606 607
25
Variable ‘importance’ 1
0.8
0.6
0.4
0.2
0
-0.1
4.5
-0.1
-0.2
0.0
3.3
0.2
2.0
0.1
0.2
3.5
-4.6
4.8
Std. effect size:
abiotic variables nitrate N
PEN (6 cm)
PEN (8 cm)
trait abund. nirK copies
biotic variables trait diversity idiosyncratic nirS evenness nirS genotype 134
R2 0.76 0.73 0.70 0.68 0.64 0.64 0.62 0.50 0.28
AICc 71.65 74.84 80.98 83.67 88.93 89.69 91.94 98.22 116.56
Δ AICc 0.00 3.19 9.33 12.02 17.28 18.04 20.29 26.57 44.91
wi 0.82 0.17 0.01 0.00 0.00 0.00 0.00 0.00 0.00
Denitrification (ng N2 O g−1 soil min−1)
●
R 2 = 0.15
20
R 2 = 0.012
20
●
●
15
●
●
● ●
●
● ●● ●
● ●
●
●
●
● ● ●
●
●● ● ●●
●
0.2
0.3
5
●
0.4
0.5
●
0.6
● ●
● ●
●
0.7
0.8
0.9
5
1.0
●
10
● ●
●
●● ●
●●
● ● ●
●
● ●
●
●
5 7.0
●
●●
●
7.2
7.4
● ●
●
● ●
● ●
● ●
●
●
7.6
●
10 ●●
● ●● ● ●
7.8
●
● ●
● ●
8.0
ni rK gene copy number g−1 soil
5
●
0.75
●● ●
●
0.80
●
●
●
●
●
●
0.7
0.85
ni rS evenness
●
● ●●
0.8
0.9
1.0
1.1
1.2
●
R 2 = 0.281
●
15 ●
●
● ●
● ● ● ● ● ● ●●●
● ●
● ● ●
●
0.90
● ● ● ●
5
●
● ● ●
● ●
●
● ● ● ●
● ● ● ●
0
●
●
●
10
● ●
●
●
●
●
● ●
●
●
● ● ● ● ●● ● ●
●
●
●
●
Penetration resistance at 8 cm (N)
●
●
●
●
●
●
15
● ●●
●
●
●
● ●
●
●
●
●
0.6
20
●
●
● ● ●
1.1
●
R 2 = 0.209
20
● ●
●
Penetration resistance at 6 cm (N)
●
15
●
●
●
●● ● ●
● ●
●
●
● ●
0.6
●
●
●
●
● ●
●
●
● ●
10
●
●
● ●●
0.5
●
R 2 = 0.087
●
●●
mg NO3 − N 100 g−1 soil 20
● ● ● ● ● ●
●
●
● ●
●
● ●
●
●
●
● ●
●
● ●
●
●
●●
● ●
●
● ●
10
●● ●●
● ●
●
● ●
●
● ● ●● ●
●
●
●
● ●
10
●
15
●
●
●
● ●●
●
●
●
15
● ●
●
R 2 = 0.003
20
●
5
Denitrification (ng N2 O g−1 soil min−1)
●
●● ● ● ● ● ● ● ● ●
● ●
●
1
2
3
4
ni rS genotype 134 percent abundance
●
2
R = 0.152
(ng N2 O g−1 soil min−1)
Denitrification
20
●
15
● ●
●
●●
10 5
●
●
●
● ●● ● ● ● ● ●
● ●●● ● ● ●● ● ● ● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ●
●
0.0
0.1
0.2
●
0.3
PEN8cm − PEN6cm (N)
0.4
