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 min class="Chemical">crobial n class="Chemical">communities are critical to carbon (C) flow in the biosphere, ecosystem models have only recently begun to simulate variation in their metabolism in soils (Li ; Wang ; Weider ). Increasingly powerful and available data on microbial community structure and function might help to better inform these efforts (McGuire and Treseder, 2010; Singh ; Schimel and Schaeffer, 2012; Graham ). Yet, compared with taxonomically constrained processes such as 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 ).

n class="Chemical">Phosphorus (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 linpan>ked to n class="Disease">biomass N: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 stoin class="Chemical">chiometry, min class="Chemical">crobial communities and soil C cycling are not currently well understood. In culture, microbial growth rates are linked with cell N:P stoichiometry and ribosomal RNA content or gene copy number (Makino and Cotner, 2004; Karpinets ; Keiblinger ; Vieira-Silva and Rocha, 2010; Franklin ), which can vary among bacterial lineages (Mouginot ; Roller ). Separately, variation or manipulation of 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 en class="Chemical">cosystem-scale relationships between soil C cycling and P availability may arise due to changes in microbial metabolism, which reflect a dependence of growth rates on 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 n class="Chemical">C flux data were obtainpan>ed from a ca. 120 hectare area under continuous rice cultivation since 2008 on Twitchell Island (38.1053 N, 121.6542 W) in the Sacramento-San Joaquin Delta, California, USA. Study plots (4 × 5 m2) were established to test N fertilization effects across a gradient of soil C content in Rindge mucky silt loam soils (Typic Halposaprists). Additional details of this study and site, including crop yield and plant N uptake, are described by Espe ). Fertilized plots had solid 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 effen class="Chemical">cts of added N onpan> tran class="Chemical">ce 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 CO2, 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 ean class="Chemical">ch plot an class="Chemical">cross the soil C gradient, CO2 fluxes were integrated over the inundated period of the 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 ean class="Chemical">ch of the experimental plots between rows of n class="Species">rice at the end of the growing season (9–10 August 2013), just prior to water drawdown and timed to coincide with peak flooded soil 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 ean class="Chemical">ch plot were n class="Chemical">composited and analyzed for pH, chloride and 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 extran class="Chemical">cted from 0.5 g of homogenized n class="Chemical">composite sample from each plot using a MoBio Laboratories (Carlsbad, CA, USA) PowerLyzer 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 aligning to large subunpan>it ribosomal RNA genes of bacterial, archaeal, fungal and other eukaryotic origins were identified and classified using the METAXA2 pipeline (Bengtsson-Palme ).

Ban class="Chemical">cterial and arn class="Chemical">chaeal 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 n class="Chemical">communpan>ity funpan>n class="Chemical">ction by PICRUSt (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

Statistin class="Chemical">cal relationpan>ships were evaluated usinpan>g 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 CO2 fluxes and C turnover were screened using linear regression, and selected relationships were evaluated using standardized major axis regression, with tests of N fertilization interactions (R: SMATR 3 package). Soil C chemistry (by Fourier-transformed infrared spectroscopy) was compared with nutrients and C turnover using partial least squares regression (R: pls package).

Relationships between n class="Chemical">C turnpan>over and the abunpan>dance of metagenome functions and taxonomic categories (aggregated at each taxonomic rank) were evaluated using 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 Unpan>iFran class="Chemical">c 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

Min class="Chemical">crobial genera n class="Chemical">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 PICRUSt-imputed pathways of C, N and P substrate metabolism (Supplementary Table S3) and ribosome copy numbers.

We n class="Chemical">calculated the relative ‘allocation’ of genes for utilizing available or recalcitrant C, N and P substrates for each OTU in our soils, based on ratios of metabolic categories (Supplementary Table S3). Each OTU was scored by its ratio of imputed copy numbers of 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 funn class="Chemical">ctionpan>al traits (inpan>cluding ribosome copy number) and microbial community composition were evaluated using a linear model for community assembly by trait selection (CATS) with the mvabund R package (Warton ). The 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 funn class="Chemical">ctionpan>al ‘guilds’ were formed usinpan>g K-means clustering (R: vegan) of preselected functional genes related to C, N and P substrate metabolism (Supplementary Table S3). Comparisons of microbial community structure (16S UniFrac distance) and de novo functional ‘guilds’ as 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 matrices. 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 signifin class="Chemical">cantly alonpan>g the soil n class="Chemical">C gradient, with effects of N fertilization appearing contingent on ambient nutrient availability. The highest cumulative CO2 fluxes occurred at intermediate levels of soil C (Figure 1a, Supplementary Figure S1), while soil C turnover (Σ g 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 n class="Chemical">classes of n class="Chemical">C, N and P elemental metabolism in shotgun 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 n class="Chemical">C turnpan>over inpan> minpan>eral soils were asson class="Chemical">ciated with more abundant carbohydrate degradation genes and inversely correlated with the abundance of aromatic compound utilization and small carboxylic (TRAP) 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).

Metabolin class="Chemical">c funpan>n class="Chemical">ctions predicted based on 16S profiles (using PICRUSt) were largely consistent (Spearman ρ2=0.82) with shotgun 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 funpan>gal abunpan>dances based on large subunit rRNA sequences in our metagenome data were generally low (0.9–3.9 and 0.1–1.4%, respectively; Supplementary Table S2).

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

n class="Species">Rice soil min class="Chemical">crobial community membership at the end of the growing season was not significantly altered by N fertilization at any level of soil C (Supplementary Figures S7a and d), from the perspective of community structure (UniFrac distances), taxonomic groups 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, pH and 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 linpan>eages were n class="Chemical">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). Phyla that were negatively correlated with soil C turnover included Euryarchaeota, Gemmatimonadetes and 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 inn class="Chemical">creasinpan>gly positive n class="Chemical">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 carbohydrate utilization and dissimilatory 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

Comparisonpan> of imputed genome allocation to C, N and P substrates (C, N and 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.

Min class="Chemical">crobial traits for n class="Chemical">C, N and P allocation and ribosome copy number each had significant interactions with soil chemical variables in 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 funn class="Chemical">ctionpan>al groups based onpan> 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 carbohydrate allocation (C), considerably improved correlations between soil C turnover and microbial ‘guilds’ (Table 4), 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, n class="Chemical">C cycling rates appeared dependent on the supply of soil C, N and P, yet constraints imposed by 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 n class="Chemical">C turnpan>over was most closely associated with P availability, including N: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, den class="Chemical">clinpan>inpan>g n class="Chemical">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 carboxylic compounds were more abundant in higher C organic soils (Figure 1f). These shifts were mirrored by changes in microbial C substrate utilization potential in our 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 effen class="Chemical">cts onpan> n class="Chemical">C turnover, N fertilization did not significantly alter microbial metabolic potential or community structure in our rice soils (Figure 2, Supplementary Figure S7a and d), in contrast to upland soils where N fertilization consistently reduces soil C turnover and alters microbial communities (Ramirez ; Fierer ; Ramirez ). Both 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).

P fertilization effen class="Chemical">cts onpan> den class="Chemical">composition also appear to be contingent on the availability of labile C (Fisk ; Nottingham ; 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

Min class="Chemical">crobial n class="Chemical">community structure varied with several factors along the soil gradient, including strong independent relationships with soil P stoichiometry (Table 1), in agreement with our second hypothesis. The abundances of several microbial groups were correlated with soil C turnover, including members of the Actinobacteria, Firmicutes, Chloroflexi, some clades of the Acidobacteria and narrower groups of 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, min class="Chemical">crobial taxa n class="Chemical">correlated with C turnover rates were enriched in imputed carbohydrate utilization genes while harboring fewer genes for 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 allon class="Chemical">cationpan> to resourn class="Chemical">ce 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 carbohydrates had less allocation to 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).

Allon class="Chemical">cationpan> to use of n class="Chemical">C, N and P substrates differed among high-level taxonomic groups (Figure 4, Supplementary Figure S8,Table 3), in agreement with our third hypothesis, that genomic features and potential utilization of C, N and 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 inpan> resource allocation strategies may underlie patterns of community assembly across environmental gradients (Weber and Deutsch, 2010; Edwards , b; van Bodegom ; Litchman , b). Our CATS models (Warton ) revealed interactions between soil nutrient availability and genomic resource allocation (Figure 4d) that suggest differences in nutrient utilization shape community composition. Genomic allocation to 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 P sn class="Chemical">cavenginpan>g allocation was less clearly related to soil nutrient availability, ribosome copy numbers had significant positive interactions with 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).

Inn class="Chemical">corporationpan> of funpan>n class="Chemical">ctional 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 P use. Accounting for differences in 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 n class="Chemical">C turnpan>over an class="Chemical">cross a large soil C:N:P gradient is in part mediated by 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 n class="Species">metagenome sequences are available from the Joinpan>t Genome Institute’s IMG database (http://img.jgi.doe.gov) unpan>der IDs 3300005873–3300005905. Shotgunpan> n 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.