A genomic perspective on stoichiometric regulation of soil carbon cycling.

Similar to plant growth, soil carbon (C) cycling is constrained by the availability of nitrogen (N) and phosphorus (P). We hypothesized that stoichiometric control over soil microbial C cycling may be shaped by functional guilds with distinct nutrient substrate preferences. Across a series of rice fields spanning 5-25% soil C (N:P from 1:12 to 1:70), C turnover was best correlated with P availability and increased with experimental N addition only in lower C (mineral) soils with N:P⩽16. Microbial community membership also varied with soil stoichiometry but not with N addition. Shotgun metagenome data revealed changes in community functions with increasing C turnover, including a shift from aromatic C to carbohydrate utilization accompanied by lower N uptake and P scavenging. Similar patterns of C, N and P acquisition, along with higher ribosomal RNA operon copy numbers, distinguished that microbial taxa positively correlated with C turnover. Considering such tradeoffs in genomic resource allocation patterns among taxa strengthened correlations between microbial community composition and C cycling, suggesting simplified guilds amenable to ecosystem modeling. Our results suggest that patterns of soil C turnover may reflect community-dependent metabolic shifts driven by resource allocation strategies, analogous to growth rate-stoichiometry coupling in animal and plant communities.

Introduction

Although microbial communities are critical to class="Chemical">carbon (C) flow in the biosclass="Chemical">phere, ecosystem models have only recently begun to simulate variation in their metabolism in soils (Li ; Wang ; Weider ). Increasingly class="Chemical">powerful and available data on microbial community structure and function might helclass="Chemical">p to better inform these efforts (McGuire and Treseder, 2010; Singh ; Schimel and Schaeffer, 2012; Graham ). Yet, comclass="Chemical">pared with taxonomically constrained class="Chemical">processes such as class="Chemical">pan class="Chemical">ammonia or methane oxidation (Bouskill ; Ho ), delineation of functional groups for decomposition in soils is challenging due to the broad distribution and functional redundancy of the relevant traits (Allison and Martiny, 2008; Schimel and Schaeffer, 2012; Berlemont and Martiny, 2013; Martiny ). However, consideration of interactions between C and nutrient cycling by soil microbes may help clarify ecologically relevant functional guilds. For example, consistent shifts in terrestrial soil C turnover and retention with mineral N addition (Craine ; Treseder, 2008) appear mediated in part by increased abundance of microbes utilizing simple C substrates, at the expense of those mining complex aromatic C and organic N (Treseder ; Fierer ; Ramirez ; Cederlund ; Amend ).

class="Chemical">Phosphorus (class="Chemical">pan class="Chemical">P) availability can also contribute to regulation of soil C cycling, and the underlying traits for microbial P cycling could analogously be linked with traits for C and N acquisition. Although less commonly studied than soil C and N interactions, P availability may affect land C sink strength at ecosystem and global scales (Wang ; Vitousek ; Goll ) and contribute to the regulation of soil C turnover rates even in nominally N-limited habitats, such as grasslands, temperate forests and leaf litter (Bradford ; Manzoni ; Strickland ; Fisk ). Soil P availability may not only affect the biomass of soil microbes (Griffiths ; Zhang ) but critically might also control community-scale rates of metabolism (Strickland ; Spohn and Chodak, 2015), which could reflect underlying stoichiometric constraints at the scale of individual cells (Hartman and Richardson, 2013).

As posited by the growth rate hypothesis (Elser , 2000, 2003), cellular growth rates are linked to class="Disease">biomass N:class="Chemical">pan class="Chemical">P ratios by the high P demands of ribosomal RNA, which determines in part the rate of synthesis of N-rich proteins. Differentiation of organism C, N and P demands based on growth rate variation forms the foundation of Ecological Stoichiometry theory (Sterner and Elser, 2002; Vrede ; Allen and Gillooly, 2009), enabling community shifts under different nutrient regimes to be connected with predictable alterations in ecosystem C cycling, particularly in aquatic ecosystems (Sterner and Elser, 2002; Weber and Deutsch, 2010; Follows and Dutkiewicz, 2011; Hessen ; Mock ). Development of a parallel framework to link stoichiometric regulation of microbial metabolism to soil C cycling at the community and ecosystem scales is highly desirable (Hall ; Sistla and Schimel, 2012; Zechmeister-Boltenstern ), especially given the integration of stoichiometric regulation of primary producers and decomposition into current terrestrial ecosystem models (Yang ; Reed ).

However, the inter-relationships between soil stoichiometry, microbial communities and soil C cycling are not currently well understood. In culture, microbial growth rates are linked with cell N:class="Chemical">P stoichiometry and ribosomal RNA content or gene coclass="Chemical">py number (Makino and Cotner, 2004; Karclass="Chemical">pinets ; Keiblinger ; Vieira-Silva and Rocha, 2010; Franklin ), which can vary among bacterial lineages (Mouginot ; Roller ). Seclass="Chemical">parately, variation or maniclass="Chemical">pulation of class="Chemical">pan class="Chemical">P availability and stoichiometry in soils has been associated with shifts in microbial community composition (Güsewell and Gessner, 2009; Fanin ; Leff ; Spohn ), without considering relationships to soil C cycling or metabolic differences among responsive microbes.

We postulated that ecosystem-scale relationships between soil C cycling and class="Chemical">P availability may arise due to changes in microbial metabolism, which reflect a declass="Chemical">pendence of growth rates on class="Chemical">pan class="Chemical">P availability, and are underpinned by shifts in microbial communities based on differences in their nutrient utilization potential. We combined plot-scale manipulations and measurements of C cycling across a broad soil C:N:P gradient with high-throughput sequencing techniques to address specific hypotheses at different scales of ecological organization, including: (1) Soil C cycling rates are associated with P availability and stoichiometry at the ecosystem scale; (2) At the community scale, taxonomy and function of microbes are associated with soil P stoichiometry; (3) At the individual scale, genomic features and potential utilization of C, N and P substrates differ characteristically among taxa, and (4) Across these scales, rates of soil C cycling are linked with the abundance of microbial groups defined by differences in stoichiometric resource allocation. To enable evaluation of the latter two hypotheses and interrogate relationships across scales, we inferred the stoichiometric utilization potential of individual microbes in situ, using 16S rRNA and shotgun metagenome sequencing in combination with inference from existing reference genomes (Langille ).

Materials and methods

Experimental design and measurement of soil C cycling and chemistry

Soils and C flux data were obtained from a ca. 120 hectare area under continuous pan class="Species">rice cultivation since 2008 on Twitchell Island (38.1053 N, 121.6542 W) in the Sacramento-San Joaquin Delta, California, USA. Study class="Chemical">plots (4 × 5 m2) were established to test N fertilization effects across a gradient of soil C content in Rindge mucky silt loam soils (Tyclass="Chemical">pic Halclass="Chemical">posaclass="Chemical">prists). Additional details of this study and site, including croclass="Chemical">p yield and class="Chemical">plant N uclass="Chemical">ptake, are described by Esclass="Chemical">pe ). Fertilized class="Chemical">plots had solid class="Chemical">pan class="Chemical">urea pellets broadcast at a single time point within 5 days of flooding on 15 May 2013.

A subset of these plots was used to test effects of added N on trace greenhouse gas fluxes in a complete randomized block design with four replicate plots each for control (0 N) and 80 kg N ha−1 fertilization (80 N) treatments in each of four fields with ca. 5, 10, 20 and 25% soil C content (32 total plots). Fluxes of class="Chemical">CO2, class="Chemical">pan class="Chemical">CH4 and N2O in closed vented soil chambers were measured weekly for the calendar year 2013, with complete methods and data described in a companion study (Ye ).

For each plot across the soil C gradient, class="Chemical">CO2 fluxes were integrated over the inundated class="Chemical">period of the class="Chemical">pan class="Species">rice-growing season (15 May–10 August; Supplementary Figure S1) to match soil collection date, instead of using closely correlated annual fluxes (Ye ). We then calculated soil C turnover, the proportion of soil C content (in the top 15 cm of soil) respired over the growing season as CO2 (Σ g CO2 m−2 per g soil C m−2) to account for differences in bulk C supply on soil respiration.

Soil samples were obtained from each of the experimental plots between rows of class="Species">rice at the end of the growing season (9–10 August 2013), just class="Chemical">prior to water drawdown and timed to coincide with class="Chemical">peak flooded soil class="Chemical">pan class="Chemical">CO2 and CH4 emissions (Supplementary Figure S1, Ye ). Three replicate topsoil samples were taken from each plot 0–15 cm beneath the soil surface, corresponding approximately to the soil horizon (anthraquic horizon) affected by agricultural activity in paddy soils (Kögel-Knabner ), which also harbors the most abundant populations of bacteria (Lee ).

Soils from each plot were composited and analyzed for pH, class="Chemical">chloride and class="Chemical">pan class="Chemical">sulfate concentrations, total pools of C, N and P and extractable NH4+, NO3− and PO42−. Further details on soil collection, processing, biogeochemistry and greenhouse gas fluxes are provided in Supplementary Methods and Supplementary Table S1. The relative abundances of C compound classes were determined in lyophilized soils using Fourier-transformed infrared spectroscopy (attenuated total reflection Fourier-transformed infrared spectroscopy), with 64 scans by a Thermo-Scientific (Waltham, MA, USA), Nicolet iS50 FT-IR and peak assignments were made following Parikh ).

Microbial community sequencing and data processing

Soil DNA was extracted from 0.5 g of homogenized composite sample from each plot using a MoBio Laboratories (Carlsbad, CA, USA) class="Chemical">PowerLyzer class="Chemical">pan class="Chemical">PowerSoil DNA Extraction Kit. For metagenomic shotgun sequencing, ~300 bp insert indexed libraries (Kapa Biosystems, Wilmington, MA, USA) were pooled 12 per lane on an Illumina (San Diego, CA, USA) HiSeq 2500 using TruSeq SBS v4 chemistry to produce 150 bp length paired-end reads, yielding ~5.0 Gbp per sample after contaminant and quality filtering (Supplementary Methods, Supplementary Table S2). Sequences were annotated via the Integrated Microbial Genomes (IMG) pipeline (Huntemann ) (details in Supplementary Methods). Gene abundance data calculated using the clusters of orthologous group (COG) ontology (Galperin ) and adjusted for contig depth were downloaded from JGI’s IMG database (Markowitz ) as a gene count by sample table from samples stored under IMG Genome IDs 3300005873–3300005905. Across our samples, ca. 42% of genes were assigned to COGs (Supplementary Table S2); others were not considered further.

Reads apan class="Chemical">ligning to large subunit ribosomal RNA genes of bacterial, archaeal, fungal and other eukaryotic origins were identified and classified using the METAXA2 class="Chemical">piclass="Chemical">peline (Bengtsson-class="Chemical">pan class="Chemical">Palme ).

Bacterial and archaeal community composition was determined from soil DNA extracts by amplifying the V4 region of 16S rRNA using barcoded multiplexed primers 515F and 816R (Caporaso ; Tremblay ). Sequences were processed using a standard pipeline (Tremblay ) yielding ca. 24 000 unique operational taxonomic units (OTUs) at depth of >74 000 reads per sample (details in Supplementary Methods).

For imputation of community function by class="Chemical">PICRUSt (class="Chemical">pan class="Chemical">Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) (Langille ), 16S rRNA sequences were mapped at 97% sequence identity to the GreenGenes reference database (DeSantis ) for which PICRUSt provides precomputed functional profiles, matching ca. 46% of sample reads and 48% of abundant OTUs (>250 total reads) (Supplementary Table S2). Abundance data from 16S sequences were corrected for ribosome copy number bias before predicting the abundance of COG ontology functions in soil metagenomes using PICRUSt. Imputed and sequenced metagenome COG function tables were rarified separately and then both normalized to 10 million COGs/sample to compare accuracy of PICRUSt predictions using the evaluate_test_datsets.py function.

Statistical data analysis

Statistical relationships were evaluated using packages in the R statistical program (R Core Team, 2015) and Quantitative Insights Into Microbial Ecology (QIIME) software (Caporaso ). Relationships between log-transformed soil biogeochemical data, soil pan class="Chemical">CO2 fluxes and C turnover were screened using linear regression, and selected relationshiclass="Chemical">ps were evaluated using standardized major axis regression, with tests of N fertilization interactions (R: SMATR 3 class="Chemical">package). Soil C chemistry (by Fourier-transformed infrared sclass="Chemical">pectroscoclass="Chemical">py) was comclass="Chemical">pared with nutrients and C turnover using class="Chemical">partial least squares regression (R: class="Chemical">pls class="Chemical">package).

Relationships between C turnover and the abundance of class="Species">metagenome functions and taxonomic categories (aggregated at each taxonomic rank) were evaluated using class="Chemical">pan class="Chemical">Pearson’s correlations with control of false discovery rates (FDRs) at ⩽5% false positive (FDR⩽0.05). Shifts in the abundances of these functions were shown in heatmaps (R:Heatmap2) using log-transformed COG count data, which were standardized (z-score scaled) across all soil plots. Aggregated shifts in C, N and P metabolism (categorical hierarchy in Supplementary Table S3) were compared with soil N:P ratios using log, Box–Cox or polynomial transformations prior to standardized major axis regression. Phylogenetic patterns were visualized using the Interactive Tree of Life web interface (Letunic and Bork, 2016).

Abundance weighted UniFrac distances (Hamady ) were used to evaluate relationships between phylogenetic community structure, soil chemistry and C turnover using Mantel’s tests, while N fertilization effects on community UniFrac distance were tested by analysis of similarities, each with 999 permutations of samples (R: vegan). Fertilization responses of log-transformed OTU and taxonomic group abundances were tested with indicator species analysis (R: vegan) and using analysis of variance models with FDR⩽0.05 (QIIME: category_significance.py).

Trait differentiation and gradient responses of imputed function

Microbial genera correlated with soil C turnover (FDR<0.05) were recoded as either positively or negatively correlated, and shrinkage linear discriminant analysis (R: sda) was used to test differences among these groups in pan class="Chemical">PICRUSt-imclass="Chemical">puted class="Chemical">pathways of C, N and class="Chemical">pan class="Chemical">P substrate metabolism (Supplementary Table S3) and ribosome copy numbers.

We calculated the relative ‘allocation’ of genes for utilizing available or recalcitrant C, N and class="Chemical">P substrates for each OTU in our soils, based on ratios of metabolic categories (Suclass="Chemical">pclass="Chemical">plementary Table S3). Each OTU was scored by its ratio of imclass="Chemical">puted coclass="Chemical">py numbers of class="Chemical">pan class="Chemical">carbohydrate to aromatic and carboxylic utilization functions (C), while nitrogen utilization was scored on assimilatory vs dissimilatory processes (N) following Helton ), and P functions were scored based on ratios of scavenging (for example, phosphatases) vs uptake (for example, P transporters) genes (P). Details of functional genes and categories can be found in Supplementary Methods and Supplementary Table S3. Patterns of C, N and P substrate allocation for each OTU were plotted on the three axes C, N and P with confidence intervals for major taxonomic groups (R: plot3Drgl). ADONIS non-parametric modeling (R: vegan) was used to test the sequential contribution of taxonomic ranks to the three-dimensional C, N and P allocation of microbial OTUs, using standardized and centered scores.

Relationships between these functional traits (including ribosome copy number) and microbial community composition were evaluated using a linear model for community assembly by trait selection (pan class="Species">CATS) with the mvabund R class="Chemical">package (Warton ). The class="Chemical">pan class="Species">CATS model represents taxon count data as a function of environmental and trait variables, with significant trait by environment interactions (filtered by LASSO) indicating ‘trait selection’.

De novo functional ‘guilds’ were formed using K-means clustering (R: vegan) of preselected functional genes related to C, N and pan class="Chemical">P substrate metabolism (Suclass="Chemical">pclass="Chemical">plementary Table S3). Comclass="Chemical">parisons of microbial community structure (16S UniFrac distance) and de novo functional ‘guilds’ as class="Chemical">predictors of trends in soil C cycling were made using Mantel’s tests to determine correlations between differences in community or guild structure and C turnover using Bray–Curtis dissimilarity matclass="Chemical">pan class="Species">rices. Guild abundances were also weighted (multiplied) by their mean ribosome copy number and/or C allocation values prior to additional Mantel’s tests to evaluate how accounting for these traits might improve relationships between guild abundances and soil C turnover.

Results

Soil C cycling mediated by N:P availability and C quality

Soil respiration varied significantly along the soil C gradient, with effects of N fertilization appearing contingent on ambient nutrient availability. The highest cumulative class="Chemical">CO2 fluxes occurred at intermediate levels of soil C (Figure 1a, Suclass="Chemical">pclass="Chemical">plementary Figure S1), while soil C turnover (Σ g class="Chemical">pan class="Chemical">CO2−C m−2 per g soil C m−2) decreased with soil C content and total N:P ratios (Figures 1b and d), consistent with global-scale patterns of higher N:P ratios with increasing soil C content (Figure 1c). Soil C turnover was also positively correlated with extractable inorganic P (Figure 1e,Table 1). N fertilization increased both CO2 production and C turnover only in lower C soils (5% and 10% C; here termed ‘mineral’), where N:P ratios were less than the ratio of 16:1 by Redfield (1958) (Figures 1a, b and d). Soil C compounds also varied with C content along the soil gradient, with mineral soils containing relatively more carbohydrates and fewer aromatic and carboxylic compounds (Figure 1f, Supplementary Figure S2) than higher C soils (20 and 25% C; here termed ‘organic’).

Figure 1

C, N and P interactions shape soil CO2 flux and soil C turnover. Effects of soil C content and added N on cumulative growing season CO2 fluxes (a) and C turnover (b) from rice soils. (c) Soil N:P ratios as a function of soil C in rice study plots (colors) and globally (gray). Soil C turnover as a function of (d) total N:P ratios and (e) soil available (Olsen) P. Correlation coefficients (R2) and dashed lines in (b, d and e) show fit of ANCOVA models (x+N treatment) using Standardized Major Axis (Type II) regression. (f) Fourier-transformed infrared spectra of soil organic matter, framed by high (dried leaf tissue) and low (clay) C reference materials, colored by field% C as in other panels. Spectra from control soils are shown as solid lines while N-fertilized treatments are shown as dashed lines. Asterisks in (a, b, d and e) show differences (analysis of variance, P⩽0.05) among control (0 N) and 80 kg N ha−1 treatments (80 N) tested separately for each rice field.

Table 1

Statistical relationships among soil C turnover, chemistry and prokaryotic community composition (by UniFrac distance)

Response	pH	Total			Extractable			Ratios
		C	N	P	NO3−	NH4+	PO42−	C:N	C:P	N:P	N:PExtr.
C turnover	0.22	0.58	0.47	—	0.29	0.38	0.62	0.65	0.70	0.69	0.72
UniFrac	0.67	0.74	0.67	0.28	0.56	0.46	0.48	0.45	0.64	0.67	0.68
partial|pH, C	—	—	—	0.32	—	—	0.43	0.40	0.41	0.41	0.43
C.turn|UniFrac	—	0.24	—	—	—	—	0.49	0.55	0.55	0.54	0.59

Data for soil C turnover responses are R2 values of log-linear standardized major axis regression models, while microbial community (UniFrac) responses are Mantel’s r correlation coefficients. Boldface numbers: correlation coefficients >0.5, and only relationships with R2 or r >0.2 and P⩽0.05 are shown. Italicized numbers indicate predictors that are not independent from C turnover by their definition, and boldface was removed.

Shifts in C, N and P substrate metabolism associated with C turnover

Major classes of C, N and class="Chemical">P elemental metabolism in shotgun class="Chemical">pan class="Species">metagenome data varied with N:P ratios along the soil C gradient (Figure 2a, Table 2), with steeper shifts in the relative abundances of aggregated functions for C, N and P metabolism appearing near the Redfield ratio (N:P=16:1). With increasing soil N:P ratios (and soil C content), dominant metabolic pathways shifted from carbohydrate to aromatic/carboxylic metabolism and from N dissimilation to assimilation, while P scavenging potential increased (phosphatases, phytases and phosphonate utilization; Figure 2a, Supplementary Figures S3–5). Conversion of these data into ratios representing preference for more or less readily available substrates (see Materials and methods section) revealed increasing use of more resistant C, N and P forms with increasing N:P along the soil C gradient (Supplementary Figures S6a–c, Supplementary Table S4). C, N and P substrate utilization patterns were also intertwined, with carbohydrate degradation linked to denitrification and P transport and aromatic/carboxylic utilization linked with N assimilation and P scavenging (Supplementary Figures S6d–f, Supplementary Table S5).

Figure 2

Microbial C, N and P substrate use linked with stoichiometry and soil C turnover. Metabolic classes share color schemes in both panels, which compare their abundances in shotgun sequence and PICRUSt imputed data for each study plot. (a) Abundance of aggregated metabolic functional classes for C (top), N (middle) and P cycling varied strongly with soil N:P ratios. Vertical scales (aggregate ortholog counts) are displayed for each metabolic category, with additional scales for PICRUSt-imputed data given in small italics where needed to align data trends but absent if data were on the same scale as metagenome sequence counts. (b) Relative abundance of specific metabolic functions correlated with soil C turnover (Pearson’s r⩾0.7; FDR<0.05) among rice soils, with both positively and negatively correlated functions shown. Additional functions for C, N and P metabolism that were not correlated with soil C are shown in Supplementary Figures S3–S5.

Table 2

Regression relationships between soil N:P ratios and aggregate categories of C, N and P substrate metabolism in metagenome and PICRUSt data (shown in Figure 2a)

Data source	Ccarb	Carom	Ndiss	Nassim	Puptk	Pscav
Metagenome	0.84	0.72	0.64	0.80	0.58	0.84
PICRUSt	0.58	0.49	0.34	0.53	0.55	0.55
Metagenome vs PICRUSt	0.70	0.55	0.48	0.34	0.17	0.81

Abbreviation: PICRUSt, Phylogenetic Investigation of Communities by Reconstruction of Unobserved States. Data are R2 from linear models, with transformation by Box–Cox or quadratic regression models (italicized) as needed to satisfy model assumptions. All models shown had P⩽0.001; boldface shows R2⩾0.5.

Higher rates of soil C turnover in mineral soils were associated with more abundant class="Chemical">carbohydrate degradation genes and inversely correlated with the abundance of aromatic comclass="Chemical">pound utilization and small carboxylic (TRAclass="Chemical">pan class="Chemical">P) transporter genes (Figure 2b, Supplementary Figure S3). Soil C turnover was also correlated with the metagenomic abundance of dissimilatory nitrate reduction and negatively correlated with P uptake and P scavenging (alkaline and acid phosphatases, phytase) genes (Figure 2b, Supplementary Figures S3–S5).

Metabolic functions predicted based on 16S profiles (using class="Chemical">PICRUSt) were largely consistent (Sclass="Chemical">pearman ρ2=0.82) with shotgun class="Chemical">pan class="Species">metagenome annotations at the plot scale (Supplementary Table S6), and similarly predicted higher C turnover was linked with more genes for carbohydrate utilization but fewer for carboxylic/aromatic utilization or P scavenging (Figure 2a), though not all patterns were reproduced, particularly for N metabolism.

Total eukaryotic and fungal abundances based on large subunit rRNA sequences in our pan class="Species">metagenome data were generally low (0.9–3.9 and 0.1–1.4%, resclass="Chemical">pectively; Suclass="Chemical">pclass="Chemical">plementary Table S2).

Microbial communities varied with C content and turnover, not N addition

pan class="Species">Rice soil microbial community membershiclass="Chemical">p at the end of the growing season was not significantly altered by N fertilization at any level of soil C (Suclass="Chemical">pclass="Chemical">plementary Figures S7a and d), from the class="Chemical">persclass="Chemical">pective of community structure (UniFrac distances), taxonomic grouclass="Chemical">ps or individual OTUs. Changes in microbial community structure were correlated with each chemical factor along the soil C gradient but were most closely aligned with soil C, class="Chemical">pH and class="Chemical">pan class="Chemical">P stoichiometry (Table 1), including soil C:P, N:P and extractable inorganic N:P. Partial Mantel’s tests excluding variance attributable to main study factors of soil C and (co-varying) pH showed significant, independent relationships between microbial community structure and soil stoichiometry, total P and inorganic P (Table 1) but not total or inorganic N forms. Soil C turnover was moderately correlated with shifts in microbial community structure (UniFrac distance; Mantel’s r=0.42), with residual variance attributable to soil stoichiometry and inorganic P (Table 1).

The abundances of several lineages were correlated with soil C turnover (Figure 3a, Supplementary Figure S7b and c, Supplementary Table S7) including Actinobacteria (primarily genera within the order Actinomycetales), Chloroflexi (Anaerolineae and Ktedonobacteria) and Fibrobacteres. The phyla Actinobacteria, Chloroflexi and Firmicutes also had significantly more OTUs that were positively correlated with C turnover than were negatively correlated (Figure 3a, Supplementary Figures S7b and c). class="Chemical">Phyla that were negatively correlated with soil C turnover included Euryarchaeota, Gemmatimonadetes and class="Chemical">pan class="Chemical">Planctomycetes (except Isosphaeraceae).

Figure 3

Imputed C, N and P substrate use and genomic features vary with phylogeny and among OTUs correlated with soil C turnover. (a) Phylogenetic patterns of soil microbes associated with soil C turnover (shown at OTU level and higher taxonomic ranks, R2⩾0.5, P⩾0.05) with imputed genome features and potential use of different C, N and P substrates. Separate color legends are used for microbial taxa, their relationships with C turnover and functional categories of C, N and P metabolism; functional categories are colored as in Figure 1. (b) Imputed ribosome copy numbers were enriched among OTUs that were more strongly correlated with C turnover. (c) Functions enriched in OTUs correlated with higher (and lower) C turnover based on linear discriminant analysis. Bars showing significantly different traits are colored to match functional categories in Figure 1 and (a).

Imputed microbial traits associated with soil C turnover

OTUs with increasingly positive correlations to soil C turnover tended to have higher imputed ribosome copy numbers (Figure 3b), and linear discriminant analysis revealed that ribosome copies were higher on average in organisms positively correlated with C turnover (Figure 3c). In parallel to community-scale functional patterns (Figure 2), microbes that were positively correlated with C turnover had more abundant class="Chemical">carbohydrate utilization and dissimilatory class="Chemical">pan class="Chemical">nitrate reduction genes (Figure 3c). High C turnover organisms also had less P scavenging and more abundant low affinity (pit) rather than high affinity (pst) phosphate transporter genes (Figure 3c).

C, N and P substrate utilization linked to community composition and C turnover

Comparison of imputed genome allocation to C, N and pan class="Chemical">P substrates (C, N and class="Chemical">pan class="Chemical">P) showed gross differences in substrate and element use strategies among OTUs (Figures 4a and b, Supplementary Figure S8), with much of the variation attributable to differences among phylum (R2=0.46) and class (R2=0.79) level taxonomic groupings (Table 3).

Figure 4

Taxonomic differentiation in C, N and P allocation contributes to microbial community assembly. (a) Taxonomic differences in stoichiometric allocation, where log-transformed axes for allocation among C, N and P substrates are derived from ratios of broad element classes and data points are individual PICRUSt-imputed genomes scaled by total COG count and colored by taxonomic groups. (b) 95% confidence intervals for the major taxonomic groups plotted as three-dimensional ellipsoids. Gray lines in (a and b) show (log) zero intercepts for each axis, representing 1:1 ratios for each allocation axis. Relationships between traits and taxa abundance in our ‘Community assembly via trait selection’ model (Warton ) illustrated in (c) were used to evaluate significant environment × trait interactions in a joint species distribution model. (d) CATS model interaction terms between stoichiometric traits and soil chemical variables for rice soils, where color indicates the sign and intensity of interactions for traits aggregated at the order level.

Table 3

Sequential contribution of taxonomic ranks to differentiation of C, N and P substrate allocation and genomic features among imputed genomes

Taxonomic rank	C:N:P allocation	Ribosome copy no.	COG count	Ccarb:arom allocation	Ndiss:assim allocation	Pscav:uptk allocation
Phylum	0.46	0.50	0.22	0.67	0.44	0.28
Class	0.79	0.68	0.49	0.87	0.80	0.71
Order	0.86	0.78	0.73	0.88	0.87	0.82
Family	0.91	0.90	0.85	0.94	0.94	0.87
Genus	0.95	0.96	0.96	0.96	0.96	0.93
OTU (resid.)	1.00	1.00	1.00	1.00	1.00	1.00

Abbreviations: COG, clusters of orthologous group; OTU, operational taxonomic unit. Data are R2 values from ADONIS models of genomic features, partitioned sequentially by taxonomic rank levels. Taxonomic contributions to C:N:P allocation were calculated from distance matrices using C, N and P substrate allocation ratios jointly. Correlation coefficients were calculated using all taxa bins at each rank, by comparing Euclidean dissimilarity matrices of standardized z-scores of log-transformed C, N and P allocation ratios.

Microbial traits for C, N and class="Chemical">P allocation and ribosome coclass="Chemical">py number each had significant interactions with soil chemical variables in class="Chemical">pan class="Species">CATS ‘trait selection’ models (Warton ) of the joint abundance of microbial orders across plots (Figures 4c and d). At each taxonomic rank, ribosome copy numbers were positively associated with inorganic P concentrations (Figure 4d, Supplementary Figure S9) and had weak negative association with nitrate. Soil nitrate concentrations had positive interactions with relative allocation to carbohydrates (C), N dissimilation (N) and, to a lesser extent, P scavenging (P). Accounting for variation in these traits improved CATS community assembly models considerably compared with null models using only chemical predictors (Supplementary Table S8).

De novo functional groups based on clustering of similar patterns of potential C and nutrient use alone were not more predictive of differences in soil C cycling rates than were differences in 16S profiles (Table 4). However, weighting the abundance of these groups by simple growth rate-related traits, including ribosome copy number and pan class="Chemical">carbohydrate allocation (C), considerably imclass="Chemical">proved correlations between soil C turnover and microbial ‘guilds’ (Table 4), class="Chemical">particularly where N was added to soils (R2=0.75).

Table 4

Comparison of de novo stoichiometric guilds with different trait weighting factors and phylogenetic community structure as predictors of soil C turnover

Predictor data	All plots	Control (no N)	Fertilized (+80 N)
CNP guilds	0.43	0.36	0.71
CNP guilds × rrn	0.46	0.38	0.74
CNP guilds × Ccarb:arom	0.48	0.44	0.74
CNP guilds × rrn ×  Ccarb:arom	0.51	0.49	0.75
16S UniFrac	0.42	0.41	0.68

Values shown are Mantel’s correlation coefficients (r) of Bray–Curtis distances among guild abundances (or UniFrac distance) and soil C turnover among plots. Stoichiometric guilds (CNP guilds) representing functionally similar organisms were delineated using K-means clustering (n=20) on imputed sets of genes for C, N and P metabolism. Guilds were weighted by their mean trait value for ribosome copy numbers (rrn) and C. The highest correlation coefficient for each set of fertilization treatment groups is highlighted in boldface.

Discussion

Stoichiometric nutrient limitation of soil C cycling and metabolism

In the wetland soils studied, C cycling rates appeared dependent on the supply of soil C, N and class="Chemical">P, yet constraints imclass="Chemical">posed by class="Chemical">pan class="Chemical">P availability featured most prominently. Non-linear patterns of soil CO2 fluxes reflected declining soil C turnover rates with increasing organic matter (Figures 1a and b). One explanation for this decline in soil C turnover is the lower availability of P relative to N as soil C accumulates (Figure 1c), consistent with global-scale patterns (Hartman and Richardson, 2013).

In agreement with our first hypothesis, soil C turnover was most closely associated with class="Chemical">P availability, including N:class="Chemical">pan class="Chemical">P ratios and inorganic P (Table 1,Figures 1d and e). Moreover, we observed a breakpoint in N fertilization effects roughly corresponding to the Redfield ratio of N:P=16:1 (Figure 1d), which delineates a shift from N to P limitation (Redfield, 1958; Sterner and Elser, 2002; Cleveland and Liptzin, 2007). Similar to interacting effects of N and P additions on soil C cycling in upland soils (Nottingham ; Poeplau ), our results suggest N fertilization altered C cycling only where N was limiting relative to P. Importantly, these results are opposite of those expected from C:N ratios alone, which would predict N fertilization effects only in high C:N organic soils.

Alternately, declining C turnover at higher soil C may reflect accumulation of aromatic compounds, which are both resistant to degradation and hypothesized to directly inhibit decomposition (Melillo ; Freeman ). Accordingly, we found aromatic and class="Chemical">carboxylic compounds were more abundant in higher C organic soils (Figure 1f). These shifts were mirrored by changes in microbial C substrate utilization class="Chemical">potential in our class="Chemical">pan class="Species">metagenome data, which revealed more carbohydrate-degrading genes in mineral soils and greater potential utilization of complex C forms in organic soils (Figure 2, Supplementary Figure S6). Notably, aromatic compounds were most closely associated with increasing C:N, C:P and N:P ratios, opposite of inorganic P availability, while increasing soil C content was linked with greater abundance of both carbohydrates and aromatics (Supplementary Figure S2).

Despite its effects on C turnover, N fertilization did not significantly alter microbial metabolic potential or community structure in our class="Species">rice soils (Figure 2, Suclass="Chemical">pclass="Chemical">plementary Figure S7a and d), in contrast to uclass="Chemical">pland soils where N fertilization consistently reduces soil C turnover and alters microbial communities (Ramirez ; Fierer ; Ramirez ). Both class="Chemical">pan class="Chemical">nitrous oxide fluxes (Ye ) and gene abundances for N dissimilation were higher in our mineral soils (Figure 2), which could suggest excess N availability consistent with their lower C:N ratios. However, N addition increased C turnover only in mineral soils (Figure 1) and did not alter in situ N2O fluxes in any of our soils (Ye ). This decoupling of denitrification from N supply could indicate that denitrifiers are instead limited by available C, as suggested by correlated N dissimilation and carbohydrate genes across the gradient (Figure 2, Supplementary Figure S6).

class="Chemical">P fertilization effects on decomclass="Chemical">position also aclass="Chemical">pclass="Chemical">pear to be contingent on the availability of labile C (Fisk ; Nottingham ; class="Chemical">pan class="Chemical">Poeplau ). The strong relationship we observed between inorganic P availability and C turnover (Figure 1) in our rice soils resembles coupling of available P with microbial metabolic rates across land use and global habitat types, including wetlands (Strickland ; Hartman and Richardson, 2013). We hypothesized that this variation in C cycling and metabolism arises from underlying shifts in microbial community structure and function, and that coupling of ecosystem biogeochemical cycles reflects metabolic coupling of element use within and among microbial taxa.

Microbial communities and traits linked with stoichiometry and C cycling

Microbial community structure varied with several factors along the soil gradient, including strong independent relationships with soil class="Chemical">P stoichiometry (Table 1), in agreement with our second hyclass="Chemical">pothesis. The abundances of several microbial grouclass="Chemical">ps were correlated with soil C turnover, including members of the Actinobacteria, Firmicutes, Chloroflexi, some clades of the Acidobacteria and narrower grouclass="Chemical">ps of class="Chemical">pan class="Chemical">Proteobacteria and Bacteroidetes (Figure 3, Supplementary Figure S7b–c, Supplementary Table S7). Many of these lineages have been previously associated with higher rates of soil C cycling or microbial metabolism (Fierer ; Goldfarb ; Fierer ; Ramirez ) potentially reflecting their greater carbohydrate utilization potential (Berlemont and Martiny, 2013; Berlemont ).

In our soils, microbial taxa correlated with C turnover rates were enriched in imputed class="Chemical">carbohydrate utilization genes while harboring fewer genes for class="Chemical">pan class="Chemical">P scavenging (Figure 3), consistent with biogeochemical and metabolic shifts observed along the gradient (Figures 1 and 2). Organisms that were highly correlated with soil C turnover also had greater imputed ribosome gene copy numbers (Figure 3), which are associated with increased growth rates among cultured microbes (Klappenbach ; Shrestha ; Vieira-Silva and Rocha, 2010; Roller ). This key trait also integrates growth rates with cellular demand for N and P (Elser , 2003), underlying functional differentiation and trophic status in higher organisms (Sterner and Elser, 2002; Wright ; Arrigo, 2005; Hessen ).

Trait differentiation underlies community assembly and C cycling

Tradeoffs in allocation to resource acquisition and growth contribute to metabolic specialization and niche partitioning among organisms (Tilman, 1990; Johnson ; Edwards ; Litchman , b; Díaz ). In our data, organisms inferred to use proportionately more class="Chemical">carbohydrates had less allocation to class="Chemical">pan class="Chemical">P scavenging (Figure 4, Supplementary Figure S8, Supplementary Table S9), consistent with their opposing patterns among OTUs positively or negatively correlated with C turnover (Figure 3), and with shifts in ecosystem C and P cycling in soil biogeochemical and metagenomic data (Figures 1 and 2, Supplementary Figure S6).

Allocation to use of C, N and class="Chemical">P substrates differed among high-level taxonomic grouclass="Chemical">ps (Figure 4, Suclass="Chemical">pclass="Chemical">plementary Figure S8,Table 3), in agreement with our third hyclass="Chemical">pothesis, that genomic features and class="Chemical">potential utilization of C, N and class="Chemical">pan class="Chemical">P substrates would differ characteristically among taxa. Taxonomic patterns in carbohydrate use and P scavenging in our data also appeared broadly similar to those observed in comparative genomic and metagenomic studies (Berlemont and Martiny, 2013; Berlemont ; Chai ). Deeply rooted partitioning of nutrient allocation could help account for consistent biogeographic patterns in microbial communities (Phillipot ; Ramirez ) and their responses to fertilization at high taxonomic ranks (Leff ; Spohn ). Higher taxonomic ranks also accounted for most of the covariation of microbial communities with chemistry in our soils (Supplementary Table S10), consistent with results compared across terrestrial and aquatic ecosystems (Lu ).

Differences in resource allocation strategies may underlie patterns of community assembly across environmental gradients (Weber and Deutsch, 2010; Edwards , b; van Bodegom ; Litchman , b). Our class="Species">CATS models (Warton ) revealed interactions between soil nutrient availability and genomic resource allocation (Figure 4d) that suggest differences in nutrient utilization shaclass="Chemical">pe community comclass="Chemical">position. Genomic allocation to class="Chemical">pan class="Chemical">carbohydrate use and N dissimilation both had positive interactions with nitrate concentrations across taxonomic ranks (Figure 4, Supplementary Figure S9), suggesting again that higher denitrification rates (Ye ) and genetic potential in our mineral soils (Figure 2) were associated with greater carbohydrate use, analogous to coupled mineral N and carbohydrate use in upland soils (Ramirez ; Amend ).

Although class="Chemical">P scavenging allocation was less clearly related to soil nutrient availability, ribosome coclass="Chemical">py numbers had significant class="Chemical">positive interactions with class="Chemical">pan class="Chemical">inorganic P concentrations at each taxonomic rank tested (Figure 4d, Supplementary Figure S9). This result could suggest that inorganic P availability limits the abundance of high ribosome copy number organisms, whose high growth rates and low C use efficiency (Roller ) may in part account for greater C turnover in our P-rich mineral soils. This mechanism, if observed elsewhere, could in part account for large-scale relationships between inorganic P availability and microbial metabolic rates across soils (Strickland ; Hartman and Richardson, 2013).

Incorporation of functional trait variation can enable community-driven metabolic modeling of nutrient-mediated changes in ecosystem C cycling (Weber and Deutsch, 2010; Bouskill ; Litchman ; Martiny ), which may be simplified by accounting for tradeoffs among traits (Edwards ; Johnson ; Litchman ). We tested simplified representations of microbial communities by deriving de novo functional ‘guilds’ of microbes based on similar patterns of imputed C, N and class="Chemical">P use. Accounting for differences in class="Chemical">pan class="Chemical">carbohydrate allocation (C) and ribosome copy numbers among these guilds yielded stronger relationships with soil C cycling in our soils than phylogenetic relationships alone (Table 4), while using fewer free parameters (Supplementary Table S11). This finding supports our fourth hypothesis that rates of soil C cycling are linked with the abundance of microbial groups defined by differences in stoichiometric resource allocation, although this awaits confirmation in other ecosystems.

Conclusions

We observed that C turnover across a large soil C:N:class="Chemical">P gradient is in class="Chemical">part mediated by class="Chemical">pan class="Chemical">P availability, reflecting differences in microbial P metabolism related to growth rates and nutrient use. Interactions among ecosystem C, N and P cycling were mirrored by community metabolic potential and appear related to shifts in community structure arising from trait-based environmental filtering of organisms by soil nutrient availability. Our findings suggest that stoichiometric resource allocation and ribosome copy number are key traits that underpin patterns of microbial community structure and function in some soils and may mediate ecosystem C cycling under multiple nutrient limitation. Grouping microbes by their resource allocation patterns improved relationships with soil nutrients and C turnover in the rice paddy soils studied, and the tradeoffs we observed among traits might help to simplify future efforts to assimilate microbial sequence data into ecosystem models.

Data availability

Annotated pan class="Species">metagenome sequences are available from the Joint Genome Institute’s IMG database (httclass="Chemical">p://img.jgi.doe.gov) under IDs 3300005873–3300005905. Shotgun class="Chemical">pan class="Species">metagenome and 16S rRNA data are deposited together in the NCBI Sequence Read Archive, under SRA study numbers SRP099852–SRP099861, SRP099869–SRP099878, SRP099883–SRP099890 and SRP099897–SRP099900.