Unearthing the Ecology of Soil Microorganisms Using a High Resolution DNA-SIP Approach to Explore Cellulose and Xylose Metabolism in Soil.

We explored microbial contributions to decomposition using a sophisticated approach to DNA Stable Isotope Probing (SIP). Our experiment evaluated the dynamics and ecological characteristics of functionally defined microbial groups that metabolize labile and structural C in soils. We added to soil a complex amendment representing plant derived organic matter substituted with either (13)C-xylose or (13)C-cellulose to represent labile and structural C pools derived from abundant components of plant biomass. We found evidence for (13)C-incorporation into DNA from (13)C-xylose and (13)C-cellulose in 49 and 63 operational taxonomic units (OTUs), respectively. The types of microorganisms that assimilated (13)C in the (13)C-xylose treatment changed over time being predominantly Firmicutes at day 1 followed by Bacteroidetes at day 3 and then Actinobacteria at day 7. These (13)C-labeling dynamics suggest labile C traveled through different trophic levels. In contrast, microorganisms generally metabolized cellulose-C after 14 days and did not change to the same extent in phylogenetic composition over time. Microorganisms that metabolized cellulose-C belonged to poorly characterized but cosmopolitan soil lineages including Verrucomicrobia, Chloroflexi, and Planctomycetes.

Introduction

Soils worldwide contain 2300 Pg of carbon (C) which accounts for nearly 80% of the C present in the terrestrial biosphere (Batjes, 1996; Amundson, 2001). Soil microorganisms drive C flux through the terrestrial biosphere and C respiration by soil microorganisms produces annually tenfold more n class="Chemical">CO2 than fossil fuel emissions (Chapin, 2002). Despite the contribution of microorganisms to global C flux, many global C models ignore the diversity of microbial physiology (Six et al., 2006; Allison et al., 2010; Treseder et al., 2011) and we still know little about the ecophysiology of soil microorganisms. Characterizing the ecophysiology of microbes that mediate C decomposition in soil has proven difficult due to their overwhelming diversity. Such knowledge should assist the development and refinement of global C models (Neff and Asner, 2001; Bradford et al., 2008; McGuire and Treseder, 2010; Wieder et al., 2013).

Though microorganisms mediate 80–90% of the soil C-cycle (Coleman and Crossley, 1996; Nannipieri et al., 2003), and microbial community composition can account for significant variation in C mineralization (Strickland et al., 2009), terrestrial C-cycle models rarely consider the community composition of soils (Zak et al., 2006; Reed and Martiny, 2007). Variation in microbial community composition can be linked effectively to rates of soil processes when diagnostic genes for specific functions are available (e.g., n class="Chemical">nitrogen fixation; Hsu and Buckley, 2009). However, the lack of diagnostic genes for describing soil-C transformations has limited progress in characterizing the contributions of individual microorganisms to decomposition. Remarkably, we still lack basic information on the physiology and ecology of the majority of organisms that live in soils. For example, contributions to soil processes remain uncharacterized for cosmopolitan bacterial phyla in soil such as Acidobacteria, Chloroflexi, Planctomycetes, and Verrucomicrobia. These phyla combined can comprise 32% of soil microbial communities (based on surveys of the SSU rRNA genes in soil; Buckley and Schmidt, 2002; Janssen, 2006).

To predict whether and how biogeochemical processes vary in response to microbial community structure, it is necessary to characterize functional niches within soil communities. Functional niches defined on the basis of microbial physiological characteristics have been successfully incorporated into biogeochemical process models (e.g., Wieder et al., 2013; Kaiser et al., 2014). In some C-cycle models physiological parameters such as growth rate and substrate specificity are used to define functional niche behavior (Wieder et al., 2013). However, it is challenging to establish the phylogenetic breadth of functional traits. Functional traits are often inferred from the distribution of diagnostic genes across genomes (Berlemont and Martiny, 2013) or from the physiology of isolates cultured on laboratory media (Martiny et al., 2013). For instance, the wide distribution of the glycolysis operon in microbial genomes is interpreted as evidence that many soil microorganisms participate in glucose turnover (McGuire and Treseder, 2010). However, the functional niche may depend less on the distribution of diagnostic genes across genomes and more on life history traits that allow organisms to compete for a given substrate as it occurs in the environment. For instance, rapid resuscitation and fast growth are traits that could allow certain microorganisms to compete effectively for n class="Chemical">glucose in environments that exhibit high temporal variability. Alternatively, metabolic efficiency and slow growth rates could be traits that allow other microbes to compete effectively for glucose in environments characterized by low temporal variability in glucose supply. These different competitive strategies would not be apparent from genome analysis, or when strains are grown in isolation. Hence, life history traits, rather than genomic capacity for a given pathway, are likely to constrain the diversity of microbes that metabolize a given C source in the soil under a given set of conditions. Therefore, to generate an understanding of functional niche as it relates to biogeochemical processes in soils it is important to characterize microbial functional traits as they occur in situ or in microcosm experiments.

Nucleic acid stable-isotope probing (SIP) links genetic identity and activity without the need of diagnostic genetic markers or cultivation and has expanded our knowledge of microbial processes (Chen and Murrell, 2010). n class="Chemical">Nucleic acid SIP has notable complications, however, including the need to add large amounts of labeled substrate (Radajewski et al., 2000), label dilution resulting in partial labeling of nucleic acids (Radajewski et al., 2000), the potential for cross-feeding and secondary label incorporation (DeRito et al., 2005a), and variation in genome G + C content (Buckley et al., 2007). As a result, most applications of SIP have targeted specialized microorganisms [for instance, methylotrophs (Lueders et al., 2004b), syntrophs (Lueders et al., 2004a), or microorganisms that target pollutants (DeRito et al., 2005b)]. Exploring the soil-C cycle with SIP has proven to be more challenging because SIP has lacked the resolution necessary to characterize the specific contributions of individual microbial groups to the decomposition of plant biomass. High throughput DNA sequencing technology, however, improves the resolving power of SIP (Aoyagi et al., 2015). It is now possible to use far less isotopically labeled substrate resulting in more environmentally realistic experimental conditions. It is also possible to sequence SSU rRNA genes from numerous density gradient fractions across multiple samples thereby increasing the resolution of a typical nucleic acid SIP experiment (Verastegui et al., 2014). With this improved resolution the activity of more soil microorganisms can be assessed. Further, since microbial activities can be more comprehensively assessed, we can begin to determine the ecological properties of functional groups defined by a specific activity in a DNA-SIP experiment. We have employed such a high resolution DNA stable isotope probing approach to explore the assimilation of both xylose and cellulose into bacterial DNA in an agricultural soil.

We added to soil a complex amendment representative of organic matter derived from fresh plant biomass. All treatments received the same amendment but the identity of isotopically labeled substrates was varied between treatments. Specifically, we set up a control treatment where all components were unlabeled, a treatment with n class="Chemical">13C-xylose instead of unlabeled xylose, and a treatment with 13C-cellulose instead of unlabeled cellulose. Soil was sampled at days 1, 3, 7, 14, and 30 and we identified microorganisms that assimilated 13C into DNA at each point in time. We designed the experiment to identify soil bacteria that metabolize easily mineralizable C (xylose) and structural C (cellulose) components of plant derived organic matter, to characterize aspects of their ecological niche, and to characterize the temporal dynamics of xylose and cellulose metabolism in soil.

Methods

Soil collection and preparation

We collected soils from an organically managed agricultural field in Penn Yan, n class="Chemical">New York. The agricultural field site has been described previously (Berthrong et al., 2013). Soils were Honoeye/Lima, a silty clay loam on calcareous bedrock, neutral pH, with 12.15 (± s.d. 0.78) mg C g−1 dry soil and 1.16 (± s.d. 0.13) mg N g−1 dry soil (Berthrong et al., 2013). On November 21, 2011, twelve soil cores (5 cm diameter × 10 cm depth) were collected in duplicate from six different sampling locations in the field using a slide hammer bulk density sampler (coordinates: (1) N 42° 40.288′ W 77° 02.438′, (2) N 42° 40.296′ W 77° 02.438′, (3) N 42° 40.309′ W 77° 02.445′, (4) N 42° 40.333′ W 77° 02.425′, (5) N 42° 40.340′ W 77° 02.420′, (6) N 42° 40.353′ W 77° 02.417′). Soils were sieved (2 mm), homogenized, distributed into flasks (10 g in each 250 ml flask, n = 34) and equilibrated at room temperature for 2 weeks until the respiration rate normalized (as determined by GCMS measurement of headspace CO2). Sieving causes a transient increase in soil respiration (Datta et al., 2014) and pre-incubation ensures that this labile organic matter is consumed and/or stabilized prior to the beginning of the experiment.

Soil microcosms

We amended soils with a complex amendment representing plant derived organic matter substituted with either n class="Chemical">13C-xylose or 13C-cellulose to represent easily mineralizable and structural C pools derived from abundant components of plant biomass. By mass the amendment contained 38% cellulose, 23% lignin (alkali extracted, Sigma), 20% xylose, 3% arabinose, 1% galactose, 1% glucose, 0.5% mannose, 10.6% amino acids (Teknova C9795), and 2.9% Murashige Skoog basal salt mixture (which contains macro and micro-nutrients associated with plant biomass; Sigma Aldrich M5524). The amendment had a C:N ratio of 10. This mixture approximates the molecular composition of switchgrass biomass with hemicellulose replaced by its constituent monomers (David and Ragauskas, 2010; Yan et al., 2010). Cellulose (2 mg cellulose g−1 d.w. soil) and lignin (1.2 mg lignin g−1 d.w. soil) were uniformly distributed over the soil surface as a powder and the remaining constituents were added in solution in a volume of 0.12 ml g−1 d.w. soil. The volume of liquid was determined in relation to soil moisture to achieve 50% water holding capacity. Water holding capacity of 50% was chosen, in relation to the texture for this soil, to achieve ~70% water filled pore space, which is the optimal water content for respiration (Linn and Doran, 1984). The total amendment was added to soil at 2.9 mg C g−1 soil dry weight (d.w.) which is representative of natural concentrations in soil during early phases of decomposition (Schneckenberger et al., 2008). The C in the amendment comprised 19 % of the total C in the soil. The cellulose-C (0.88 mg C g−1 soil d.w.) and xylose-C (0.42 mg C g−1 soil d.w.) in the amendment comprised 6 and 3% of the total C in the soil, respectively.

We set up three parallel treatments varying the isotopically labeled component in each treatment (see Figure S1 for outline of experimental design). The treatments were: (1) 12C Control Plant Simulant (12CCPS) which included only unlabeled components, (2) 13C-Cellulose Plant Simulant (n class="Chemical">13CCPS) which included 13C-cellulose (96 atom % 13C, synthesized as described below) instead of unlabeled cellulose, and (3) 13C-Xylose Plant Simulant (13CXPS) which included 13C5-D-xylose (98 atom % 13C, Sigma Aldrich) instead of unlabeled xylose. A total of 12 microcosms were established for the 13C-xylose and control treatments and 10 for the 13C-cellulose treatment. Microcosms were sampled destructively at days 1 (control and xylose only), 3, 7, 14, and 30 and soils were stored at −80°C until nucleic acid extraction. Isotopic analysis was performed at the Cornell University Stable Isotope Laboratory. While microcosm experiments are useful for evaluating microbial activity under defined experimental conditions, results from such experiments need to be verified in the field before drawing conclusions about microbial activity in the natural environment.

Cellulose production

Bacterial cellulose was produced by Gluconoacetobacter n class="Species">xylinus grown in Heo and Son (2002) minimal media (HS medium) made with 0.1% glucose and without inositol. For the production of 13C-cellulose, 13C6-D-glucose, 99 atom % 13C (Cambridge Isotope) was used. Cellulose was produced in 1L Erlenmeyer flasks containing 100 mL HS medium inoculated with three colonies of G. xylinus grown on HS agar plates. Flasks were incubated statically in the dark at 30°C for 2–3 weeks. Cellulose pellicules were decanted, rinsed with deionized water, suspended in two volumes of 1% Alconox, and then autoclaved. Cellulose pellicules were purified by dialysis for 12 h in 1 L deionized water and dialysis was repeated 10 times. Harvested pellicules were dried overnight (60°C), cut into pieces, and ground using a 5100 Mixer/Mill (SPEX SamplePrep, Metuchen, NJ), and dry sieved to 53–250 μm. The particulate size range was selected to be representative of particulate organic matter in soils.

The purity of ground cellulose was checked by biological assay, Benedict's reducing n class="Chemical">sugars assay, Bradford assay, and isotopic analysis. E. coli is not able to use cellulose as a C source but is capable of growth on a variety of nutrients available in the HS medium. The biological assay consisted of E. coli inoculated into minimal M9 media which lacked a C source and was supplemented with either: (1) 0.01% glucose, (2) 2.5 mg purified, ground cellulose, (3) 25 mg purified, ground cellulose, (4) 25 mg purified, ground cellulose and 0.01% glucose. Growth in media was checked by spectrometer (OD450). No measurable growth was observed with either 2 mg or 25 mg cellulose, indicating absence of contaminating nutrients that can support growth of E. coli. In addition, the presence of 25 mg cellulose did not inhibit the growth of E. coli cultures provided with glucose (relative to control), indicating the absence of compounds in the purified cellulose that could inhibit microbial growth (data not shown). Purified cellulose was also assayed for residual proteins and sugars using Bradford and Benedict's assays, respectively. Bradford assay was performed as in Bradford (1976). Ground, purified cellulose contained 6.92 μg protein mg cellulose−1(i.e., 99.31% purity). Reducing sugars were not detected in cellulose using Benedict's reducing sugar assay (Benedict, 1909) tested at 10 mg cellulose ml−1 (data not shown). Finally, 13C-cellulose had an average 96% ± 5 (s.d.) degree of 13C labeling as determined by isotopic analysis (UCDavis Stable Isotope Facility).

Nucleic acid extraction

Nucleic acids were extracted from 0.25 g soil using a modified Griffiths protocol (Griffiths et al., 2000). Cells were lysed by bead beating for 1 min at 5.5 m s−1 in 2 mL lysis tubes containing 0.5 g of 0.1 mm diameter n class="Chemical">silica/zirconia beads (treated at 300°C for 4 h to remove RNAses), 0.5 mL extraction buffer (240 mM Phosphate buffer 0.5% N-lauryl sarcosine), and 0.5 mL phenol-chloroform-isoamyl alcohol (25:24:1) for 1 min at 5.5 m s−1. After lysis, 85 μL 5 M NaCl and 60 μL 10% hexadecyltriammonium bromide (CTAB)/0.7 M NaCl were added to lysis tube, vortexed, chilled for 1 min on ice, and centrifuged at 16,000 x g for 5 min at 4°C. The aqueous layer was transferred to a new tube and reserved on ice. To increase DNA recovery, the pellet was back extracted with 85 μL 5 M NaCl and 0.5 mL extraction buffer. The aqueous extract was washed with 0.5 mL chloroform:isoamyl alcohol (24:1). Nucleic acids were precipitated by addition of 2 volumes polyethylene glycol solution (30% PEG 8000, 1.6 M NaCl) on ice for 2 h, followed by centrifugation at 16,000 x g, 4°C for 30 min. The supernatant was discarded and pellets were washed with 1 mL ice cold 70% EtOH. Pellets were air dried, resuspended in 50 μL TE and stored at −20°C. To prepare nucleic acid extracts for isopycnic centrifugation as previously described Buckley et al. (2007), DNA was size selected (>4 kb) using 1% low melt agarose gel and β-agarase I enzyme extraction per manufacturers protocol (New England Biolab, M0392S). Final resuspension of DNA pellet was in 50 μL TE.

Isopycnic centrifugation and fractionation

We fractionated DNA on density gradients for n class="Chemical">13C-xylose treatments (days 1, 3, 7, 14, 30), 13C-cellulose treatments (days 3, 7, 14, 30), and control treatments (days 1, 3, 7, 14, 30). A total of 5 μg DNA was added to each 4.7 mL CsCl density gradient. Density gradient were composed of 1.69 g ml−1 CsCl ml−1 in gradient buffer solution (pH 8.0 15 mM Tris-HCl, 15 mM EDTA, 15 mM KCl). Centrifugation was performed at 55,000 rpm 20°C for 66 h using a TLA-110 rotor in a Beckman Coulter Optima MAX-E ultracentrifuge. Fractions of ~100 μL were collected from below by displacing the DNA-CsCl-gradient buffer solution in the centrifugation tube with water using a syringe pump at a flow rate of 3.3 μL s−1 (Manefield et al., 2002). Fractions were collected in Acroprep 96 filter plates (part no. 5035, Pall Life Sciences). The refractive index of each fraction was measured using a Reichart AR200 digital refractometer modified as previously described to measure a volume of 5 μL (Buckley et al., 2007). Buoyant density was calculated from the refractive index as previously described Buckley et al. (2007) using the equation ρ = aη − b, where ρ is the density of the CsCl (g ml−1), η is the measured refractive index, and a and b are coefficient values of 10.9276 and 13.593, respectively, for CsCl at 20°C (Birnie, 1978). The refractive index (Ri) was corrected to account for the Ri of the gradient buffer using the equation: Ri = Ri − (Ri − 1.3333). A total of 35 fractions were collected from each gradient and the average density between fractions was 0.0040 g mL−1. The DNA was desalted by washing with TE (5X 200 μL) in the Acroprep filter wells. DNA was resuspeneded in 50 μL TE.

PCR amplification

SSU rRNA genes were amplified from gradient fractions (n = 20 per gradient) and from non-fractionated Dn class="Chemical">NA from soil. Barcoded primers consisted of: 454-specific adapter B, a 10 bp barcode, a 2 bp linker (5′-CA-3′), and 806R primer for reverse primer (BA806R); and 454-specific adapter A, a 2 bp linker (5′-TC-3′), and 515F primer for forward primer (BA515F). Each PCR contained 1.25 U μl−1 AmpliTaq Gold (Life Technologies, Grand Island, NY; N8080243), 1X Buffer II (100 mM Tris-HCl, 500 mM KCl, pH 8.3), 2.5 mM MgCl2, 200 μM of each dNTP, 0.5 mg ml−1 BSA, 0.2 μM BA515F, 0.2 μM BA806R, and 10 μL of 1:30 DNA template in 25 μl total volume. The PCR conditions were 95°C for 5min followed by 22 cycles of 95°C for 10 s, 53°C for 30 s, and 72°C for 30 s, followed by a final elongation at 72°C for 5 min. Amplification products were checked by 1% agarose gel. Reactions were performed in triplicate and pooled. Amplified DNA was gel purified (1% low melt agarose) using the Wizard SV gel and PCR clean-up system (Promega, Madison, WI; A9281) per manufacturer's protocol. Samples were normalized by SequalPrep normalization plates (Invitrogen, Carlsbad, CA; A10510) or based on PicoGreen DNA quantification and pooled in equimolar concentration. Amplicons were sequenced on Roche 454 FLX system using titanium chemistry at Selah Genomics (Columbia, SC).

DNA sequence analysis

SSU rRNA gene sequences were initially screened by maximum expected errors at a specific read length threshold (Edgar, 2013). Reads that had more than 0.5 expected errors at a length of 250 nt were discarded. The remaining reads were aligned to the n class="Disease">Silva Reference Alignment as provided in the Mothur software package using the Mothur NAST aligner (DeSantis et al., 2006; Schloss et al., 2009). Reads that did not align to the expected region of the SSU rRNA gene were discarded. After expected error and alignment based quality control. The remaining quality controlled reads were annotated using the “UClust” taxonomic annotation framework in QIIME (Caporaso et al., 2010; Edgar, 2010). We used 97% sequence identity cluster seeds from the Silva SSU rRNA database (release 111Ref; Quast et al., 2013) as reference for taxonomic annotation [provided on the QIIME website, (Quast et al., 2013)]. Quality control screening filtered out 344,472 of 1,720,480 sequencing reads leaving 1,376,008 reads for downstream analyses. Reads annotated as “Chlorloplast,” “Eukaryota,” “Archaea,” “Unassigned” or “mitochondria” were culled from the dataset.

Sequences were distributed into OTUs with a centroid based clustering algorithm [i.e., UPARSE; Edgar (2013)]. The centroid selection also included chimera screening (Edgar, 2013). OTU centroids were established at a threshold of 97% sequence identity and non-centroid sequences were mapped back to centroids. Reads that could not be mapped to an OTU centroid at greater than or equal to 97% sequence identity were discarded. OTU centroids were compared (BLAST; Altschul et al., 1990; Camacho et al., 2009) to sequences in “The All-Species Living Tree” project (LTP). The LTP (version 115) is a collection of SSU rRNA gene sequences for named and cultured species of Archaea and Bacteria (Yarza et al., 2008).

We used SSU-Align (Nawrocki et al., 2009; n class="Chemical">Nawrocki and Eddy, 2013) to align SSU rRNA gene sequences. Columns in the alignment that were aligned with poor confidence (<95% of characters had posterior probability >95%) were not considered when building the phylogenetic tree leaving a multiple sequence alignment of 216 columns. Additionally, the alignment was trimmed to coordinates such that all sequences in the alignment began and ended at the same positions. FastTree was used with default parameters to build the phylogeny (Price et al., 2010).

Statistical analyses

NMDS ordination was performed on weighted Unifrac (Lozupone and Knight, 2005) distances between samples. The n class="Chemical">Phyloseq (McMurdie and Holmes, 2013) wrapper for Vegan (Oksanen et al., 2015; both R packages) was used to compute sample values along NMDS axes. The “Adonis” function in Vegan was used to perform permutational MANOVA (default parameters; Anderson, 2001).

We used DESeq2 (R package), an RNA-Seq differential expression statistical framework (Love et al., 2014), to identify OTUs that were enriched in high density gradient fractions from n class="Chemical">13C-treatments relative to corresponding gradient fractions from control treatments [for review of RNA-Seq differential expression statistics applied to microbiome OTU count data see (McMurdie and Holmes, 2014)]. Briefly, DESeq2 includes several features that enable robust estimates of standard error in addition to reliable ranking of logarithmic fold change (LFC; i.e., gamma-Poisson regression coefficients) in OTU relative abundance even with low count OTUs. Figures S2, S3 demonstrate raw data for one example responder and non-responder OTU, respectively. Responders increase in relative abundance in the high density fractions relative to control due to 13C-labeling of their DNA. As our data is compositional, unlabeled OTUs often had consistent relative abundance across the density gradients indicating the OTU DNA concentration across the gradient mirrored that of the total DNA concentration. For each OTU, we calculated LFC and corresponding standard errors for enrichment in high density gradient fractions of 13C treatments relative to control. We define “high density gradient fractions” as gradient fractions whose density falls between 1.7125 and 1.755 g ml−1. Subsequently, a one-sided Wald test was used to statistically assess LFC values. The user-defined null hypothesis was that LFC was less than one standard deviation greater than the mean of all LFC values. We independently filtered OTUs prior to multiple comparison corrections on the basis of sparsity eliminating OTUs that failed to appear in at least 45% of high density fractions for a given comparison. P-values were adjusted for multiple comparisons using the Benjamini and Hochberg method (Benjamini and Hochberg, 1995). We selected a false discovery rate of 10% to denote statistical significance.

We identified OTUs that changed in relative abundance over time again using gamma-poisson regression (R package DESeq2, Love et al., 2014). Specifically, we used day as an ordered factor as the regressor with LFC of the relative abundance in non-fractionated DNA as the outcome in the general linear model. We used the default DESeq2 base mean independent filtering and disabled the Cook's cutoff outlier detection. The null model was that abundance did not change with time and we assessed significance at a false discovery rate of 10%.

We estimated the rrn copy number for each OTU as described; i.e., we used the code and reference information provided by the authors (Kembel et al., 2012) directly. In brief, OTU centroid sequences were inserted into a reference SSU rRNA gene phylogeny (Matsen et al., 2010) from organisms of known rrn copy number. The rrn copy number was then inferred from the phylogenetic placement in the reference phylogeny.

The Net Relatedness Index (n class="Chemical">NRI) and Nearest Taxon Index (NTI) were calculated using the “picante” R package (Kembel et al., 2010). We used the “independentswap” null model for phylogenetic distribution. Briefly, the NRI and NTI evaluate phylogenetic clustering against a null model for the distribution of a trait in a phylogeny (Webb, 2000). The NRI and NTI values are z-scores and thus the greater the magnitude of the NRI/NTI, the stronger the evidence for clustering (positive values) or overdispersion (negative values). NRI assesses basal clustering whereas the NTI assesses terminal clustering (Evans and Wallenstein, 2014). The consenTRAIT clade depth for xylose and cellulose responders was calculated using R code from the original publication describing the metric (Martiny et al., 2013) which employs the R “adephylo” package (Jombart and Dray, 2010). The consenTRAIT metric is a measure of the average clade depth for a trait in a phylogenetic tree.

We calculated the change in center of mass for buoyant density, , for the relative abundance of each OTU as the change in an OTU's density profile center of mass (see below) between corresponding control and labeled gradients (Figure S2). DNA buoyant density (BD) increases with atom % n class="Chemical">13C. Therefore, the magnitude of for an OTU in response to 13C labeling provides a measure of the degree of isotopic labeling for an OTU. Because all gradients did not span the same density range and gradient fractions cannot be taken at specific density positions, we limited our analysis to the density range for which fractions were taken for all gradients. Within this density range we linearly interpolated 20 evenly spaced relative abundance values. The center of mass for an OTU along the density gradient was then the density weighted average where weights were the linearly interpolated relative abundance values. The values determined in this way do not represent the isotopic incorporation of any specific strain but rather represent an average for each OTU. This average can be shifted by the ratio of active and inactive strains within each OTU as well as by isotope dilution effects. Furthermore, values based on compositional data can be distorted in comparison to values based on absolute DNA concentration. Hence, values calculated in this way provide a qualitative estimate of the extent of labeling (for example, does DNA labeled in response to 13C-xylose additions contain on average more 13C than DNA labeled in response to the addition 13C-cellulose) but do not provide precise quantification of isotope incorporation. Note, that these values were not used to identify OTUs that incorporated 13C into DNA (for the identification of 13C-responsive OTUs, see instead LFC calculations as described above).

All code to take raw SSU rRNA gene sequencing reads to final publication figures and through all presented analyses is located at the following URL: https://github.com/chuckpr/CSIP_succession_data_analysis. Dn class="Chemical">NA sequences are available at the NCBI Short Read Archive (Accession SRP072785).

Results

We sequenced SSU rRNA gene amplicons from a total of 277 n class="Chemical">CsCl gradient fractions representing 14 CsCl gradients and 12 non-fractionated DNA samples. The SSU rRNA gene data set contained 1,102,685 total sequences. The average number of sequences per sample was 3816 (s.d. 3629) and 265 samples had over 1000 sequences. We sequenced SSU rRNA gene amplicons from an average of 20 fractions per CsCl gradient (s.d. 0.57). The average density between fractions was 0.0040 g mL−1. The sequencing effort recovered a total of 5940 OTUs. 2943 of the total 5940 OTUs were observed in non-fractionated samples. We observed 33 unique phylum and 340 unique genus annotations. Analysis of SSU rRNA gene amplicons from non-fractionated soil DNA (Figure S4) indicated that variance in community composition was mainly driven by time and that, in the absence of SIP, little variance in community composition could be attributed to isotopic treatment or block effects (i.e., bottle to bottle variance was small relative to variance over time).

Community-level signal of 13C-assimilation in relation to substrate and time

The soil microbial community respired 65% of the added xylose within 1 day and 29% of the n class="Chemical">13C from xylose remained in the soil at day 30 (Figure 1). In contrast, cellulose was respired at an average rate of 18 μg C d−1 g−1 soil d.w. and 40% of added 13C from cellulose remained in the soil at day 30 (Figure 1).

Figure 1

Distinct dynamics of metabolization are observed for . Symbols indicate the percentage of added 13C that remains in soil over time. Soils were pooled (three samples per time point per treatment) prior to measuring 13C-content. The bars indicate enumerate OTUs found by SIP to assimilate 13C into DNA (i.e., counts of xylose and cellulose responders) over time (days 1, 3, 7, 14, 30 for xylose, and days 3, 7, 14, 30 for cellulose). Note that 13C lost from soil must be lost in gaseous form, and that 13C that remains in soil may be either unmetabolized, transformed into a new chemical form, or assimilated into microbial biomass.

Distinct dynamics of metabolization are observed for . Symbols indicate the percentage of added 13C that remains in soil over time. Soils were pooled (three samples per time point per treatment) prior to measuring n class="Chemical">13C-content. The bars indicate enumerate OTUs found by SIP to assimilate 13C into DNA (i.e., counts of xylose and cellulose responders) over time (days 1, 3, 7, 14, 30 for xylose, and days 3, 7, 14, 30 for cellulose). Note that 13C lost from soil must be lost in gaseous form, and that 13C that remains in soil may be either unmetabolized, transformed into a new chemical form, or assimilated into microbial biomass.

We initially assessed assimilation of 13C into microbial Dn class="Chemical">NA by comparing the SSU rRNA gene sequence composition of SIP density gradient fractions between 13C treatments and the unlabeled control. In the unlabeled control treatment, fraction density represented the majority of the variance in SSU rRNA gene composition (Figure 2). This result is expected because genome G + C content correlates positively with DNA buoyant density and influences SSU rRNA gene composition in gradient fractions (Buckley et al., 2007). In contrast, isotope assimilation into DNA will cause variation in gene sequence composition between corresponding density fractions from controls and labeled treatments. For example, the SSU rRNA gene composition in gradient fractions from the 13C-cellulose treatment deviated from corresponding control fractions on days 14 and 30 and this difference was observed only in the high density fractions (>1.7125 g mL−1, Figure 2). Likewise, SSU rRNA gene composition in gradient fractions from the 13C-xylose treatment also deviated from corresponding control fractions but on days 1, 3, and 7 as opposed to 14 and 30 (Figure 2). The 13C-cellulose and 13C-xylose treatments also differed from each other in corresponding high density gradient fractions indicating that different microorganisms were labeled across time in these treatments (Figure 2).

Figure 2

NMDS ordination of SSU rRNA gene sequence composition in gradient fractions shows that variation between fractions is correlated with fraction density, isotopic labeling, and time. Dissimilarity in SSU rRNA gene sequence composition was quantified using the weighted UniFrac metric. SSU rRNA gene sequences were surveyed in twenty gradient fractions at each sampling point for each treatment (Figure S1). 13C-labeling of DNA is apparent because the SSU rRNA gene sequence composition of gradient fractions from 13C and control treatments differ at high density. Each point on the NMDS plot represents one gradient fraction. The size of each point is positively correlated with density and colors indicate the treatment (A) or day (B).

NMDS ordination of SSU rRn class="Chemical">NA gene sequence composition in gradient fractions shows that variation between fractions is correlated with fraction density, isotopic labeling, and time. Dissimilarity in SSU rRNA gene sequence composition was quantified using the weighted UniFrac metric. SSU rRNA gene sequences were surveyed in twenty gradient fractions at each sampling point for each treatment (Figure S1). 13C-labeling of DNA is apparent because the SSU rRNA gene sequence composition of gradient fractions from 13C and control treatments differ at high density. Each point on the NMDS plot represents one gradient fraction. The size of each point is positively correlated with density and colors indicate the treatment (A) or day (B).

We observed further differences in the pattern of isotope incorporation over time for each treatment. For example the SSU rRNA gene sequence composition in the n class="Chemical">13C-cellulose treatment was similar on days 14 and 30 in corresponding high density fractions indicating similar patterns of isotope incorporation into DNA on these days (Figure 2). In contrast, in the 13C-xylose treatment, the SSU rRNA gene composition varied between days 1, 3, and 7 in corresponding high density fractions indicating different patterns of isotope incorporation into DNA on these days (Figure 2). In the 13C-xylose treatment on days 14 and 30 the SSU gene composition was similar to control on days 14 and 30 for corresponding high density fractions (Figure 2) indicating that 13C was no longer detectable in bacterial DNA on these days for this treatment. These results show that the dynamics of isotope incorporation into DNA varied considerably for organisms that assimilated C from either xylose or cellulose.

Temporal dynamics of OTU relative abundance in non-fractionated DNA from soil

We monitored the soil microbial community over the course of the experiment by surveying SSU rRNA genes in non-fractionated Dn class="Chemical">NA from the soil. The SSU rRNA gene composition of the non-fractionated DNA changed with time (Figure S4, P = 0.023, R2 = 0.63, permutational MANOVA with model formula:weighted unifrac distances ~ day as an ordered factor). In contrast, the microbial community could not be shown to change with treatment (P-value 0.23, permutational MANOVA with model formula:weighted unfrac distance ~ treatment; Figure S4). The latter result demonstrates that substitution of 13C-labeled substrates for unlabeled equivalents could not be shown to alter the soil microbial community composition. Twenty-nine OTUs exhibited sufficient statistical evidence (adjusted P < 0.10, Wald test) to conclude they changed in relative abundance in the non-fractionated DNA over time (Figure S5). When SSU rRNA gene abundances were combined at the taxonomic rank of “class,” the classes that changed in abundance (adjusted P < 0.10, Wald test) were the Bacilli (decreased), Flavobacteria (decreased), Gammaproteobacteria (decreased), and Herpetosiphonales (increased) (Figure S6). Of the 29 OTUs that changed in relative abundance over time, 14 putatively incorporated 13C into DNA (see below and Figure S5). OTUs that likely assimilated 13C from 13C-cellulose tended to increase in relative abundance with time whereas OTUs that assimilated 13C from 13C-xylose tended to decrease (Figure S5). OTUs that responded to both substrates did not exhibit a consistent relative abundance response over time (Figure S5).

Changes in the phylogenetic composition of 13C-labeled OTUs with substrate and time

If an OTU exhibited strong evidence for assimilating 13C into Dn class="Chemical">NA, we refer to that OTU as a “responder” (see Methods for our operational definition of “responder”). The SSU rRNA gene sequences produced in this study were binned into 5940 OTUs and we assessed evidence of 13C-labeling from both 13C-cellulose and 13C-xylose for each OTU. Forty-one OTUs responded to 13C-xylose, 55 OTUs responded to 13C-cellulose, and 8 OTUs responded to both xylose and cellulose (Figures 3, 4, Figure S7, Tables S1, S2). The number of xylose responders peaked at days 1 and 3 and declined with time (Figure 1). In contrast, the number of cellulose responders increased with time peaking at days 14 and 30 (Figure 1).

Figure 3

Enrichment of OTUs in either . Each point indicates the LFC for a single OTU. High enrichment values indicate an OTU is likely 13C-labeled. Different colors represent different phyla and different panels represent different days. The final column shows the frequency distribution of LFC values in each row. Within each panel, shaded areas are used to indicate one standard deviation (dark shading) or two standard deviations (light shading) about the mean of all LFC values.

Figure 4

Phylogenetic position of cellulose responders and xylose responders in the context of all OTUs that passed sparsity independent filtering criteria (see Methods). Only those phyla that contain responders are shown. Colored dots are used to identify xylose responders (green) and cellulose responders (blue). The heatmaps indicate enrichment in high density fractions relative to control (represented as LFC) for each OTU in response to both 13C-cellulose (13CCPS, leftmost heatmap) and 13C-xylose (13CXPS, rightmost heatmap) with values for different days in each heatmap column. High enrichment values (represented as LFC) provide evidence of 13C-labeled DNA.

Phylogenetic position of cellulose responders and n class="Chemical">xylose responders in the context of all OTUs that passed sparsity independent filtering criteria (see Methods). Only those phyla that contain responders are shown. Colored dots are used to identify xylose responders (green) and cellulose responders (blue). The heatmaps indicate enrichment in high density fractions relative to control (represented as LFC) for each OTU in response to both 13C-cellulose (13CCPS, leftmost heatmap) and 13C-xylose (13CXPS, rightmost heatmap) with values for different days in each heatmap column. High enrichment values (represented as LFC) provide evidence of 13C-labeled DNA.

The phylogenetic composition of xylose responders changed with time (Figures 3, 5) and most n class="Chemical">xylose responders shared >97% SSU rRNA gene sequence identity with bacteria cultured in isolation (Table S1). On day 1, Bacilli OTUs represented 84% of xylose responders (Figure 5) and the majority of these OTUs were closely related to cultured representatives of the genus Paenibacillus (Table S1, Figure 4). For example, “OTU.57” had a strong signal of 13C-labeling at day 1 coinciding with its maximum relative abundance in non-fractionated DNA. The relative abundance of “OTU.57” declined until day 14 and it did not appear to be 13C-labeled after day 1 (Figure 6). On day 3, Bacteroidetes OTUs comprised 63% of xylose responders (Figure 5) and these OTUs were closely related to cultured representatives of the Flavobacteriales and Sphingobacteriales (Table S1, Figure 4). For example, “OTU.14,” annotated as a flavobacterium, had a strong signal for 13C-labeling at days 1 and 3 coinciding with its maximum relative abundance in non-fractionated DNA. Its relative abundance then declined until day 14 and did not show evidence of 13C-labeling beyond day 3 (Figure 6). Finally, on day 7, Actinobacteria OTUs represented 53% of the xylose responders (Figure 5) and these OTUs were closely related to cultured representatives of Micrococcales (Table S1, Figure 4). For example, “OTU.4,” annotated as Agromyces, had signal for 13C-labeling on days 1, 3 and 7 with the strongest evidence of 13C-labeling at day 7 and did not appear 13C-labeled at days 14 and 30. The relative abundance of “OTU.4” in non-fractionated DNA increased until day 3 and then declined until day 30 (Figure 6). Proteobacteria were also common among xylose responders at day 7 where they comprised 40% of xylose responder OTUs. Notably, Proteobacteria represented the majority (6 of 8) of OTUs that responded to both cellulose and xylose (Figure S7).

Figure 5

Xylose reponders in the . The left column shows counts of 13C-xylose responders in the Actinobacteria, Bacteroidetes, Firmicutes and Proteobacteria at days 1, 3, 7, and 30. The right panel shows OTU enrichment in high density gradient fractions (gray points, expressed as fold change) for responders as well as a boxplot for the distribution of fold change values [The box extends one interquartile range, whiskers extend 1.5 times the IR, and small dots are outliers (i.e., beyond 1.5 times the IR)]. Each day in the right column shows all responders (i.e., OTUs that responded to xylose at any point in time). High enrichment values indicates OTU DNA is likely 13C-labeled.

Figure 6

Raw data from individual responders highlighted in the main text (see Results). The left column shows OTU relative abundance in density gradient fractions for the indicated treatment pair at each sampling point. Time is indicated by the line color (see legend). Gradient profiles are shaded to represent the different treatments where orange represents “control,” blue “13C-cellulose,” and green “13C-xylose.” The right column shows the relative abundance of each OTU in non-fractionated DNA. Enrichment in the high density fractions of 13C-treatments indicates an OTU likely has 13C-labeled DNA.

Xylose reponders in the . The left column shows counts of n class="Chemical">13C-xylose responders in the Actinobacteria, Bacteroidetes, Firmicutes and Proteobacteria at days 1, 3, 7, and 30. The right panel shows OTU enrichment in high density gradient fractions (gray points, expressed as fold change) for responders as well as a boxplot for the distribution of fold change values [The box extends one interquartile range, whiskers extend 1.5 times the IR, and small dots are outliers (i.e., beyond 1.5 times the IR)]. Each day in the right column shows all responders (i.e., OTUs that responded to xylose at any point in time). High enrichment values indicates OTU DNA is likely 13C-labeled.

Raw data from individual responders highlighted in the main text (see Results). The left column shows OTU relative abundance in density gradient fractions for the indicated treatment pair at each sampling point. Time is indicated by the line color (see legend). Gradient profiles are shaded to represent the different treatments where orange represents “control,” blue “13C-cellulose,” and green “n class="Chemical">13C-xylose.” The right column shows the relative abundance of each OTU in non-fractionated DNA. Enrichment in the high density fractions of 13C-treatments indicates an OTU likely has 13C-labeled DNA.

The phylogenetic composition of cellulose responders did not change with time to the same extent as the n class="Chemical">xylose responders. Both the relative abundance and the number of cellulose responders increased gradually over time peaking at days 14 and 30 (Figures 1, 3, 7). Cellulose responders belonged to the Proteobacteria (46%), Verrucomicrobia (16%), Planctomycetes (16%), Chloroflexi (8%), Bacteroidetes (8%), Actinobacteria (3%), and Melainabacteria (1 OTU; Table S2). Most cellulose responders increased in non-fractionated DNA relative abundance with time and had the strongest evidence for 13C-labeling at days 14 and 30 (e.g., “OTU.6,” “OTU.64,” “OTU.2192,” “OTU.573”; Figure 6). Also, in contrast to xylose responders, many cellulose responders were only distantly related to cultured isolates (Figure 4). For example, most of the Verrucomicrobia cellulose responders fell within unidentified Spartobacteria clades that had < 85% SSU rRNA identity to cultivated isolates (Figure 4). Likewise, Chloroflexi cellulose responders belonged to an unidentified clade within the Herpetosiphonales (Figure 4) and they shared < 89% SSU rRNA identity to cultivated isolates (Figure 4). Bacteroidetes included responders to both cellulose and xylose, though the former were members of Cytophagales while the latter belonged either to Flavobacteriales or Sphingobacteriales (Figure 4). Proteobacterial cellulose responders belonged to Alpha (13 OTUs), Beta (4 OTUs), Gamma (5 OTUs), and Delta-proteobacteria (6 OTUs). Most proteobacterial cellulose responders were closely related (>97% identity) to isolated strains, including representatives of the genera: Cellvibrio, Devosia, Rhizobium, and Sorangium, which are all known for their ability to degrade cellulose (Table S2).

Figure 7

Change in relative abundance in non-fractionated DNA over time for xylose responders (13CXPS) and cellulose responders (13CCPS). Each panel represents a responders to the indicated substrate [i.e., cellulose (13CCPS) or xylose (13CXPS)] within the indicated phylum except for the lower right panel which shows all reponders to both xylose and celluose. The abbreviations Proteo., Verruco., and Plancto., correspond to Proteobacteria, Verrucomicrobia, and Planctomycetes, respectively.

Change in relative abundance in non-fractionated DNA over time for n class="Chemical">xylose responders (13CXPS) and cellulose responders (13CCPS). Each panel represents a responders to the indicated substrate [i.e., cellulose (13CCPS) or xylose (13CXPS)] within the indicated phylum except for the lower right panel which shows all reponders to both xylose and celluose. The abbreviations Proteo., Verruco., and Plancto., correspond to Proteobacteria, Verrucomicrobia, and Planctomycetes, respectively.

Characteristics of cellulose and xylose responders

Cellulose responders, relative to n class="Chemical">xylose responders, tended to have lower relative abundance in non-fractionated DNA, demonstrated signal consistent with higher atom % 13C in labeled DNA, and had lower estimated rrn copy number (Figure 8). OTUs that assimilated C from either cellulose or xylose were also clustered phylogenetically (see below) indicating that these traits were not dispersed randomly across bacterial species.

Figure 8

Characteristics of xylose responders (green) and cellulose responders (blue) based on estimated . The estimated rrn copy number of all responders is shown vs. time (A). Kernel density histogram of values shows cellulose responders had higher average than xylose responders indicating higher average atom % 13C in OTU DNA (B). The final panel indicates the rank relative abundance of all OTUs observed in the non-fractionated DNA (C) where rank was determined at day 1 (bold line) and relative abundance for each OTU is indicated for all days by colored lines (see legend). Xylose responders (green ticks) have higher relative abundance in non-fractionated DNA than cellulose responders (blue ticks). All ticks are based on day 1 relative abundance.

Characteristics of xylose responders (green) and n class="Chemical">cellulose responders (blue) based on estimated . The estimated rrn copy number of all responders is shown vs. time (A). Kernel density histogram of values shows cellulose responders had higher average than xylose responders indicating higher average atom % 13C in OTU DNA (B). The final panel indicates the rank relative abundance of all OTUs observed in the non-fractionated DNA (C) where rank was determined at day 1 (bold line) and relative abundance for each OTU is indicated for all days by colored lines (see legend). Xylose responders (green ticks) have higher relative abundance in non-fractionated DNA than cellulose responders (blue ticks). All ticks are based on day 1 relative abundance.

In the non-fractionated DNA, n class="Chemical">xylose responders [3.5 × 10−3 (s.d. 5.2 × 10−3)] had higher relative abundance than cellulose responders [1.2 × 10−3 (s.d. 3.8 × 10−3)] (Figure 8, P = 1.12 × 10−5, Wilcoxon Rank Sum test). Consistent with this observation, seven of the ten most abundant responders were xylose responders, and six of the ten most abundant OTUs overall (i.e., responders + non-responders) were xylose responders.

DNA buoyant density (BD) increases in proportion to atom % n class="Chemical">13C. Hence, the extent of 13C incorporation into DNA can be evaluated by the difference in BD between 13C-labeled and unlabeled DNA. We calculated for each OTU its mean BD weighted by relative abundance to determine its “center of mass” within a given density gradient. We then quantified for each OTU the difference in center of mass between control gradients and gradients from 13C-xylose or 13C-cellulose treatments (see methods for the detailed calculation; Figure S2). We refer to the change in center of mass position for an OTU in response to 13C-labeling as . Cellulose responder [0.0163 g mL−1 (s.d. 0.0094)] was greater than that of xylose responders (0.0097 g mL−1 (s.d. 0.0094); Figure 8, P = 1.8610 × 10−6, Wilcoxon Rank Sum test).

We predicted the rrn gene copy number for responders as described (Kembel et al., 2012). The ability to proliferate after rapid nutrient influx correlates positively to a microorganism's rrn copy number (Klappenbach et al., 2000). Cellulose responders possessed fewer estimated rrn copy numbers [2.7 (1.2 s.d.)] than n class="Chemical">xylose responders [6.2 (3.4 s.d.)] (P = 1.878 × 10−9, Wilcoxon Rank Sum test, Figure 8, Figure S8). Furthermore, the estimated rrn gene copy number for xylose responders was inversely related to the day of first response (P = 2.02 × 10−15, Wilcoxon Rank Sum test, Figure 8, Figure S8).

We assessed the phylogenetic clustering of 13C-responsive OTUs with n class="Chemical">NTI and NRI and quantified the average clade depth of cellulose and xylose responders with the consenTRAIT metric as described in methods. NRI values indicate that cellulose responders clustered overall and at the tips of the phylogeny (NRI: 4.49, NTI: 1.43) while xylose responders clustered terminally (NRI: −1.33, NTI: 2.69). The consenTRAIT clade depth for xylose and cellulose responders was 0.012 and 0.028 SSU rRNA gene sequence dissimilarity, respectively. As reference, the average clade depth as inferred from genomic analyses or growth in culture is ~0.017, 0.013, and 0.034 SSU rRNA gene sequence dissimilarity for arabinase (arabinose like xylose is a five C sugar found in hemicellulose), glucosidase and cellulase, respectively (Berlemont and Martiny, 2013; Martiny et al., 2013). These results indicate xylose responders form terminal clusters dispersed throughout the phylogeny while cellulose responders form deep clades of terminally clustered OTUs.

Discussion

We highlight two key results with implications for understanding structure-function relationships in soils, and for applying DNA-SIP in future studies of the soil-C cycle. First, most n class="Chemical">cellulose responders belonged to undescribed taxonomic groups. This suggests we still have much to learn about the diversity of cellulolytic bacteria in soil. Second, the xylose response was characterized by a succession in activity from Paenibacillus OTUs (day 1) to Bacteroidetes (day 3) and finally Micrococcales (day 7). Notably, Paenibacillus have been previously shown by DNA-SIP to metabolize glucose (Verastegui et al., 2014). This activity succession was mirrored by relative abundance profiles and may indicate trophic-C exchange between these groups. While trophic exchange has been observed previously in DNA-SIP studies (Lueders et al., 2004b) most applications of DNA-SIP focus on proximal use of labeled substrates. However, with increased sensitivity, DNA-SIP is well-suited to tracking C flows throughout microbial communities over time. Though bacteria are often considered as a single trophic level (Moore et al., 1988), trophic interactions will critically influence how the global soil-C reservoir will respond to climate change (Crowther et al., 2015). Additionally, our results show that DNA-SIP results can change dramatically over time suggesting that multiple time points are necessary to rigorously and comprehensively describe which microorganisms consume 13C-labeled substrates in nucleic acid SIP incubations.

Bacteria that consumed 13C-cellulose were members of groups found widely in soil including Spartobacteria, Planctomycetes, and Chloroflexi, for which few related isolates are available in cultivation. Often n class="Chemical">cellulose responders were less than 90% related to their closest cultured relatives showing that we can infer little of their physiology from culture-based studies. Verrucomicrobia represented 16% of the cellulose responders. Verrucomicrobia are cosmopolitan soil microorganisms that can make up to 23% of SSU rRNA gene sequences in soils (Bergmann et al., 2011) and 9.8% of soil SSU rRNA (Buckley and Schmidt, 2001). Genomic analyses and laboratory experiments show that various isolates within the Verrucomicrobia are capable of methanotrophy, diazotrophy, and cellulose degradation (Wertz et al., 2011; Otsuka et al., 2012). Moreover, Verrucomicrobia have been hypothesized to degrade polysaccharides in many environments (Chin et al., 1999; Fierer et al., 2013; Herlemann et al., 2013). The majority of verrucomicrobial cellulose responders belonged to two clades that fell within the Spartobacteria (Figure 4). Spartobacteria outnumbered all other Verrucomicrobia phylotypes in SSU rRNA gene surveys of 181 globally distributed soil samples (Bergmann et al., 2011). Given their ubiquity and abundance in soil as well as their demonstrated incorporation of 13C from 13C-cellulose, Verrucomicrobia lineages, particularly Spartobacteria, may be important contributors to global cellulose turnover.

Other notable cellulose responders include OTUs in the Planctomycetes and Chloroflexi both of which have previously been shown to assimilate n class="Chemical">13C from 13C-cellulose added to soil (Schellenberger et al., 2010). Planctomycetes are common in soil (Janssen, 2006), comprising 4 to 7% of bacterial cells in many soils (Zarda et al., 1997; Chatzinotas et al., 1998) and 7 ± 5% of SSU rRNA (Buckley and Schmidt, 2003). Although soil Planctomycetes are widespread, their activities in soil remain uncharacterized. Plantomycetes represented 16% of cellulose responders and shared <92% SSU rRNA gene sequence identity to their most closely related cultured isolates. Chloroflexi are known for metabolically diverse lifestyles ranging from anoxygenic phototrophy to organohalide respiration (Hug et al., 2013) and are among the six most abundant bacterial phyla in soil (Janssen, 2006). Recent studies have focused on Chloroflexi roles in C-cycling (Goldfarb et al., 2011; Cole et al., 2013; Hug et al., 2013) and several Chloroflexi isolates use cellulose (Goldfarb et al., 2011; Cole et al., 2013; Hug et al., 2013). Four of the five Chloroflexi cellulose responders belong to a single clade within the Herpetosiphonales (Figure 4).

Finally, a single cellulose responder belonged to the Melainabacteria phylum (95% shared SSU rRn class="Chemical">NA gene sequence identity with Vampirovibrio chlorellavorus). The phylogenetic position of Melainabacteria is debated but Melainabacteria have been proposed to be a non-phototrophic sister phylum to Cyanobacteria. An analysis of a Melainabacteria genome (Rienzi et al., 2013) suggests the genomic capacity to degrade polysaccharides though Vampirovibrio chlorellavorus is an obligate predator of green alga (Gromov and Mamkaeva, 1972).

The turnover of cellulose and plant-derived n class="Chemical">sugars in soil has been studied previously using DNA-SIP (e.g., Verastegui et al., 2014). Similar to our study, phylotypes among the Chloroflexi, Bacteroidetes, and Planctomycetes have all been previously implicated in soil cellulose degradation (Schellenberger et al., 2010). Additionally, functional metagenomics enabled by DNA-SIP has identified glycoside hydrolases putatively belonging to Cellvibrio and Spartobacteria further suggesting a role for these organisms in cellulose breakdown in soils (Verastegui et al., 2014). Fungi undoubtedly also contribute to the decomposition of cellulose in soils (Boer et al., 2005), but they are not a focus of this experiment.

In addition to taxonomic identity, we quantified four ecological properties of microorganisms that were actively engaged in labile and structural C decomposition in our experiment: (1) time of activity, (2) estimated rrn gene copy number, (3) phylogenetic clustering, and (4) density shift in response to 13C-labeling. n class="Chemical">Xylose was consumed before cellulose and these substrates were consumed by different microorganisms (Figures 1–3). Consumers of xylose had higher estimated rrn gene copy number than cellulose consumers (Figure 8). rrn copy number is positively correlated with the ability to resuscitate quickly in response to nutrient influx (Klappenbach et al., 2000) which may be the advantage that enabled xylose responders to rapidly consume xylose. Both xylose and cellulose responders were terminally clustered in phylogenetic trees suggesting that the ability to use these substrates was phylogenetically constrained. Although labile C consumption is generally considered to be mediated by a diverse set of microorganisms, we found that xylose responders at day 1 were mainly members of one genus, Paenibacillus (Figure 4). Our results suggest that life-history traits such as the ability to resuscitate quickly and/or grow rapidly may be more important in determining the diversity of microorganisms that actually mediate a given process than the genomic potential for substrate utilization. And last, xylose consumers, in contrast to cellulose consumers, had lower in response to 13C-labeling. This result suggests that many xylose C consumers were generalists, assimilating C from a variety of sources both labeled and unlabeled, while structural C consumers were more likely to be specialists and more closely associated with C from a single source.

Responders do not necessarily assimilate 13C directly from n class="Chemical">13C-xylose or 13C-cellulose, and in many ways, knowledge of secondary C-transformations and/or microbial biomass turnover may be more interesting with respect to the soil C-cycle than knowledge of primary degradation. In the case of cellulose degradation, non-cellulolytic bacteria may gain access to 13C from 13C-cellulose via ecological interactions including mutualism, parasitism, predation, and ‘cheating’. Cheaters in particular, who benefit from the activities of cellulolytic extracellular enzymes without bearing the costs of producing these enzymes, have been hypothesized to impact rates of cellulose degradation (Allison, 2005). Cheaters should be able to use oligosaccharides and sugars produced by extracellular enzymes. We note with interest that few OTUs responded both to 13C-cellulose and 13C-xylose. This means that if these responders include taxa labeled by secondary C-transformations, then these secondary C-transformations are substrate specific.

The response to 13C -xylose suggests that n class="Chemical">xylose-C moved through different trophic levels within the soil bacterial food web. Trophic exchange would explain well the precipitous drop in abundance of Paenibacillus after day 1 with subsequent 13C-labeling of Bacteroidetes at day 3 as well as the precipitous drop in abundance of Bacteroidetes at day 3 followed by 13C-labeling of Micrococcales at day 7. The Bacilli degraded xylose first (65% of the xylose-C had been respired by day 1) representing 84% of day 1 xylose responders (Figure 5). Bacilli also comprised about 6% of SSU rRNA genes present in non-fractionated DNA on day 1 (Figure 7). However, few Bacilli remained 13C-labeled by day 3 and their abundance declined reaching about 2% of soil SSU rRNA genes by day 30 (Figures 5, 7). Members of the Bacillus (Cleveland et al., 2007) and Paenibacillus in particular (Verastegui et al., 2014) have been previously implicated as labile C decomposers. Bacteroidetes OTUs appeared 13C-labeled at day 3 concomitant with the decline in relative abundance and loss of 13C-label for Bacilli. Finally, Actinobacteria appeared 13C-labeled at day 7 as Bacteroidetes xylose responders declined in relative abundance and became unlabeled. Hence, it seems reasonable to propose that Bacteroidetes and Actinobacteria xylose responders became labeled via the consumption of 13C derived from 13C-labeled microbial biomass as opposed to primary degradation of 13C-xylose. Trophic exchange could be enabled by mother cell lysis (in the case of spore formers such as Paenibacillus), viral lysis, and/or predation.

The inferred physiology of Actinobacteria and Bacteroidetes xylose responders provides further evidence for C transfer by saprotrophy and/or predation. Certain Bacteroidetes have been shown to become n class="Chemical">13C-labeled after the addition of live 13C-labeled Escherichia coli to soil (Lueders et al., 2006) indicating their ability to assimilate C from microbial biomass. In addition, the dominant OTU labeled in the 13C-xylose treatment from the Micrococcales shares 100% SSU rRNA gene sequence identity to Agromyces ramosus a known predator that feeds upon on many microorganisms including yeast and Micrococcus luteus (Casida, 1983). Agromyces are abundant microorganisms in many soils and Agromyces ramosus was the most abundant xylose responder in our experiment—and the fourth most abundant OTU in our dataset. It is notable however, that Agromyces ramosus remains unlabeled in response to 13C-cellulose. This result suggests that if Agromyces ramosus is acting as a predator in our experiment, it is either a selective predator or it is only active under a restricted set of certain environmental conditions.

If trophic transfer caused the dynamics we observed, at least three different ecological groups exchanged C in 7 days. Models of the soil C cycle often exclude trophic interactions between soil bacteria [e.g., Moore et al. (1988)], yet when soil C models do account for predators and/or saprophytes, trophic interactions are predicted to have significant effects on the fate of soil C (Kaiser et al., 2014). Climate change is expected to diminish bottom-up controls on microbial growth increasing the importance on top-down biological interactions for mitigating positive climate change feedbacks (Crowther et al., 2015). Currently the extent of bacterial predatory activity in soil, and its consequences for the soil C-cycle and carbon use efficiency is largely unknown. Elucidating the identities of bacterial predators in soil will assist in assessing the implications of climate change on global soil-C storage.

Microorganisms govern C-transformations in soil and thereby influence global climate but we still do not know the specific identities of microorganisms that carry out critical C transformations. In this experiment, cellulose responders included members of the Verrucomicrobia (Spartobacteria), Chloroflexi, Bacteroidetes, and Planctomycetes. Spartobacteria in particular are globally cosmopolitan soil microorganisms and are often the most abundant Verrucomicrobia in soil (Bergmann et al., 2011). Our results suggest three major conclusions: (i) as yet uncharacterized cellulolytic bacteria appear to be present in several major lineages of uncultivated bacteria found widely in soil, (ii) most of the OTUs that responded to n class="Chemical">xylose were labeled as a result of trophic transfer and not primary assimilation, and (iii) life history traits are likely to act as a filter constraining the diversity of active microorganisms relative to those with the genomic potential for a given metabolism.

Author contributions

CP conceived, developed and implemented data analysis and interpreted data, AC designed and performed the experiments, and both CP and AC performed data analysis, generated figures, and wrote and edited the manuscript. CP and AC contributed equally to the manuscript and should be considered co-first authors. CK and SB assisted with sampling and experimental analyses including development and implementation of the DNA sequencing pipeline. DB conceived the experimental design and supervised all aspects of experimentation, data analysis, and manuscript preparation.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.