Soil samples were collected from three sites (Blodgett, Brandy City, and Lowell Hill) in the Sierra Nevada of California which were harvested and reforested with ponderosa pine 16 years previously as part of the Cohasset soil series of the Long-Term Soil Productivity Study (Powers, 2006). The distance between sites ranged from 20 to 57 km and all shared similar forest cover and soil type (Mesic Ultic Haploxeralfs). Four treatments were sampled at each site: an unharvested reference plot (REF) and three harvested treatments accompanied by varying degrees of OM removal. Harvested treatments consisted of OM1, where tree stems (trunks) were removed, but branches and woody debris were retained; OM2, where whole tree biomass was removed; and OM3, where whole tree biomass plus the upper organic layer of the soil were removed. Photographs of harvested treatments are displayed in Figure S1. Triplicate samples were taken at each plot (45 m²) each comprised of five sub-sampled points along a plot transect to account for heterogeneity and ensure sufficient soil material. In sampling, the litter layer was removed from vegetation-free soil and the organic layer (the O-horizon) was collected with a trowel. Next, the top 20 cm of mineral soil (including the A and occasionally upper B-horizon) was collected using a Stoney auger (5 cm diameter). Samples were stored at −80°C and processed within 1 year. For an overview of sample collection and the experimental design consult Figure 1A.

Organic and mineral layer soil samples from REF, OM1 and OM3 were incubated in microcosms with no additional substrate or with one of three milled lignocellulosic substrates derived from Douglas-fir: (i) “lignocellulose,” from debarked, untreated Douglas-fir woodchips, (ii) “lignin + cellulose,” from steam treated woodchips, where hemicellulose was solubilized and removed, and (iii) “cellulose,” from steam treated woodchips which were subsequently delignified (Kumar et al., 2012). Microcosms were prepared by adding 4.5 g dry wt soil to 30-mL serum vials, adjusting moisture content to 60% (mineral) and 125% (organic) (w/v) and pre-incubated at 20°C for 1 week. Substrate was then added (10% w/w) along with CO₂ traps, consisting of sterile glass vials containing 2 mL NaOH (1M). Microcosms were then incubated at 20°C for 14 days based on time course experiments described in Wilhelm et al. (2014). Net respiration was determined by titration of the NaOH traps according to methods described by Haney et al. (2008). OM2 samples were not tested here, and in small number of other experiments, due to limitations in the quantity of available substrate.

Soil from all samples (n = 72) was incubated in paired treatments: one amended with 10% w/w of ¹³C-labeled cellulose (99 atom % ¹³C) and another with the same amount of unlabeled cellulose (natural abundance ¹³C: ~1%). The ¹²C-control incubations were included to correct for natural ¹³C content in SIP-phospholipid fatty acid (PLFA) work and to control for native populations with higher GC content (i.e., slightly heavier DNA) in SIP-DNA work, as described below. Bacterial cellulose (>99% glucose) was produced by feeding Gluconacetobacter xylinus with ¹³C-labeled glucose (see Supplementary Data 5, for details). Microcosm preparation was identical to previously described respiration assays, except for the following differences: 1.5 g (organic) and 2 g (mineral) dry wt soil were used, and incubations were for 11 days (organic) and 14 days (mineral). Following incubation, soil was lyophilized and stored at −80°C until processing. All SIP-PLFA, SIP-pyrotag and SIP-metagenomic data were derived from the same set of microcosms. Libraries termed “in situ” were derived from corresponding field soil samples that were not incubated and, post-hoc, from publicly available pyrotag libraries from LTSP sites in British Columbia (Hartmann et al., 2012).

PLFAs were extracted according to Bligh and Dyer (1959) and the ¹³C-content was analyzed using IRMS (University of British Columbia Stable Isotope Facility) ported with gas chromatography as detailed in Churchland et al. (2013). Peak identification was based on retention time compared against two reference standards: the bacterial acid methyl-ester standard (47080-0; Sigma–Aldrich, St. Louis) and a 37-Component fatty acid methyl-ester mix (47885-U; Sigma–Aldrich, St. Louis). Unidentifiable ¹³C-enriched peaks were also included in analysis if they met the following conditions: (i) detection in 3 or more samples, (ii) average δ ¹³C > +50%0 and (iii) confirmed as long-chain alkane methyl esters by GC-MS. Taxonomic affiliations of specific PLFAs were assigned according to Högberg et al. (2013), with c18:1ω9 and c18:3ω6 added as additional fungal PLFAs (Ruess and Chamberlain, 2010).

DNA was extracted from soil (0.5 g) with the manufacturer's recommended protocol for the FastDNA™ Spin Kit for Soil (MPBio, Santa Ana, CA). The mass and the atom % ¹³C of DNA extracts were measured with UHPLC-MS/MS according to Wilhelm et al. (2014). DNA extracts from replicates within each site were pooled in equal amounts and unlabeled controls were processed identically (n = 24 × 2). ¹³C-enriched DNA was recovered by density gradient ultracentrifugation according to methods in Neufeld et al. (2007) and Wilhelm et al. (2014), with improvements for greater recovery of DNA (see Supplementary Data 5). Both SIP and in situ pyrotag libraries were prepared from the 16S rRNA gene (V1–V3 regions) as well as fungal internal transcribed spacer region (ITS2) according to the procedure of Hartmann et al. (2012). Metagenomic libraries were prepared from 40–50 ng of enriched DNA using the Nextera DNA Sample Preparation Kit (Illumina Inc., CA, USA). Four shotgun metagenome libraries were generated, ¹³C-libraries from REF, OM1 and OM3 treatments as well as a ¹²C-library from the REF treatment, by pooling the corresponding DNA extracts from mineral layer samples at all three sites. These four libraries were multiplexed on two lanes of Illumina HiSeq (2 × 100-bp), yielding 285 million paired-end reads. There was insufficient ¹³C-enriched DNA to generate metagenomes for the organic layer samples.

Statistical analyses were performed using the R platform (v. 3.1.0; R Core Team, 2015). 16S rRNA gene libraries were quality-filtered and processed using Mothur (Schloss et al., 2009) according to the Schloss “454 SOP” (http://www.mothur.org/wiki/454_SOP; accessed May 2013) and were clustered into operational taxonomic units (OTUs) at 0.01% dissimilarity. ITS libraries were processed according to Hartmann et al. (2012) to create OTUs, but, due to the hypervariability of the ITS region, were also grouped based on taxonomic classification using UNITE (Kõljalg et al., 2013). We used three methods to identify OTUs differentially abundant between ¹²C-control and ¹³C-enriched library samples: “DESeq” (Anders and Huber, 2010), “limma-voom” (Ritchie et al., 2015) and uncorrected, averaged relative abundance. An OTU was deemed ¹³C-enriched if it had at least a 3-fold higher relative abundance in ¹³C vs. ¹²C libraries according to one or more of the methods. The identification of carbohydrate-active enzyme (CAZy') genes was based on BLASTX searches using methods in Cardenas et al. (2015). The following glycosyl hydrolase (GH) families contain enzymes with endoglucanase activity: GH5, 6, 7, 8, 9, 12, 16, 44, 45, 48, 51, 61, 74, 81, and 131. Shotgun metagenome libraries were preprocessed using Trimmomatic (Bolger et al., 2014; v. 0.32), to trim sequencing primers and low quality ends, and the FastX Toolkit (Gordon and Hannon, 2010; v. 0.7), to filter short or low quality reads. Paired-end and orphaned reads were all assembled using the default setting of Ray-meta (kmer size = 39 bp) (Boisvert et al., 2012; v. 2.3.1). Subsequent binning of contigs into putative genome bins was performed with Metawatt (Strous et al., 2012; v. 2.1), based on tetranucleotide frequency, and MetaBAT (Kang et al., 2014; v. 0.18.6), based on both tetranucleotide frequency and covariance in read abundance mapped to the super assembly. The completeness of genome bins was assessed by scanning for essential single-copy, house-keeping genes with hidden Markov models provided by Albertsen et al. (2013). Taxonomic designations were based on lowest-common ancestor analysis via MEGAN using matches to the NCBI “nr” database (v. 5.10.1; Huson et al., 2007). Reads from metagenomes were mapped back to genome bins using Bowtie2 (Langmead and Salzberg, 2012) to estimate their relative abundances among harvested treatments. Additional details can be found in Supplementary Data 5, while a script for all R analyses and raw data can be found in Supplementary Data 2.

Raw sequence data were deposited at the European Nucleotide Archive under the study accession (PRJEB9761) for 16S rRNA gene pyrotags (ERS803692-ERS803739) ITS pyrotags (ERS803740-ERS803786), binned genomes (ERZ288956 - ERZ288966), and metagenomic libraries (ERS1099581- ERS1099584). Raw data used in all other analysis, namely soil chemistry data, net respiration, PLFA, CAZyme abundances (among others) can be found in Supplementary Data 2.

Comparisons of decomposer activity were based on net respiration and total ¹³C assimilated from cellulose into PLFA and DNA. Net respiration was significantly lower in mineral soil microcosms from harvested plots, OM1 [t_{(2, 103)} = −3.18, p < 0.001] and OM3 [t_{(2, 103)} = −4.99, p < 0.001], relative to reference plots (REF), with or without amendment of Douglas-fir lignocellulosic substrates (Figure 1C). OM1 soil had the highest carbon content (Table 1), consistent with the retention of woody debris, yet produced less CO₂ than REF even though carbon content was weakly correlated with respiration (Spearman's r = 0.19; p = 0.049). Respiration in cellulose-amended soil from OM3 was particularly low, suggesting that cellulose-degrading populations were disproportionately affected by the greatest degree of OM removal [t_{(11, 94)} = −2.65; p < 0.01; arrow in Figure 1C]. Respiration was also positively correlated with pH (Spearman's r = 0.34; p < 0.001), which was lowest in OM1, but not with total ¹²C PLFA (Spearman's r = 0.05; p = 0.82), total nitrogen (Spearman's r = 0.18; p = 0.07) or C:N ratio (Spearman's r = −0.07; p = 0.45). Organic layer soils respired 3-fold more CO₂ than mineral soils, but no significant differences among OM removal treatments were observed. Across all experiments, measurements of organic layer soils were highly variable. In respiration assays, differences in respiration even among organic soil microcosms with or without added substrate were obscured by this variability.

The incorporation of ¹³C into PLFAs and DNA was generally lower in soils from harvested plots, with the clearest evidence occurring in mineral soils (Table 1). The total assimilation by bacteria and fungi also differed among harvested treatments, with significantly greater bacterial assimilation in both REF and OM1 (Figure 2). The differences in bacterial and fungal activity was evident in both ¹²C and ¹³C PLFA profiles, driven by an increase in fungal biomass in OM3 (1.8-fold higher than REF) relative to Gram-positive and Gram-negative bacterial biomass in REF, 1.4-fold and 1.3-fold higher than OM3, respectively. The proportion of ¹³C-enriched Gram-positive PLFAs was highest in REF (47%) followed by OM1 (35%), OM2 (31%), and OM3 (28%). Overall, the organic layer contained 4-fold greater microbial biomass, based on total PLFAs, and exhibited greater total cellulolytic activity. Minor differences in cellulolytic activity were observed among sites, which corresponded to differences in total biomass (Table S1).

The successful targeting of cellulolytic populations in sequencing libraries via SIP was supported by the following: (i) significantly higher quantities of ¹³C in soil DNA extracts and a corresponding 2.5-fold higher concentration of DNA in heavy CsCl gradient fractions (Figure S2); (ii) distinct clustering of samples according to ¹³C-enrichment and soil layer in NMDS ordinations (Figure 3); (iii) the vastly improved assembly of metagenomes from ¹³C-enriched DNA (~20%) relative to the control library (<1%) (Table S2) and (iv) lower alpha-diversity of all ¹³C-libraries relative to ¹²C- and in situ libraries (Figure S3). Harvested treatments accounted for ~9% of the total variation in ¹³C-pyrotag OTU profiles (perMANOVA; F = 1.3; p = 0.04), which was comparable to the amount explained by soil layer (Figure S4). Harvesting was not a significant factor in explaining differences among fungal ¹³C-pyrotag profiles. Harvesting did not produce differences in alpha diversity (Shannon-Wiener diversity; Figure S3) or beta diversity (UniFrac) of ¹³C-pyrotag libraries, suggesting no substantial loss or gain of cellulolytic groups occurred. However, harvesting did alter the relative abundance of major taxa incorporating ¹³C from cellulose during incubations.

¹³C-pyrotag libraries from harvested plots had diminished relative abundances of putatively cellulolytic Verrucomicrobia (Chthoniobacter and unclassified Opitutaceae), uncl. Streptomycetaceae, Burkholderia, uncl. Rhizobiaceae and Caulobacter (Table 2). Populations of Verrucomicrobia, Streptomycetaceae (Kitasatospora sp.) and Caulobacteraceae were sufficiently abundant and active (i.e., differentially abundant in ¹³C-libraries) to recover sizeable metagenome-assembled genomes (MAGs) (Table S3). Further evidence for the contraction of these populations as a result of harvesting was found in the reduced proportion of metagenomic reads which mapped to their MAGs (Figure 4) as well as their reduced relative abundances in ¹²C- and in situ pyrotag libraries (Figure 5). Conversely, a number of putatively cellulolytic taxa increased in abundance in soils from harvested plots, including a number of Betaproteobacteria and members of Myxococcales, Planctomycetes and fungi belonging to the ascomycotal family Chaetomiaceae and basidiomycotal genus Clitopilus (Table 2). The active member of Chaetomiacae was also sufficiently abundant and active to recover its 43-Mb MAG (Table S3). Overall, Ascomycota were between ~10³ and 10⁴ times more abundant than Basidiomycota in ¹³C-pyrotag libraries. The ratio of Basidiomycota to Ascomycota did not significantly differ among harvested plots in either ¹²C- or ¹³C-libraries; however, the ratio did significantly decrease in situ in OM2 [t_{(3, 65)} = −2.97; p < 0.01] and OM3 [t_{(3, 65)} = −4.18; p < 0.01] relative to REF (Figure 6).

The increased ratio of fungi to Gram-positive bacteria in harvested treatments, observed in PLFA data, was corroborated by similar changes in the relative abundance of Streptomycetaceae (Actinobacteria) and Chaetomiaceae (Ascomycota) in pyrotag libraries and shotgun metagenomes. The trend was apparent in both ¹³C- and in situ pyrotag libraries (Figure 5) as well as by read mapping to the Kitasatospora (Streptomycetaceae) and Myceliophthora thermophila MAGs (Chaetomiaceae) (Figure 4). The decline of Streptomycetaceae (Kitasatospora), along with Caulobacter and Opitutaceae, and increase in relative abundance of Chaetomiaceae was corroborated by previously published pyrotag libraries from LTSP field sites in British Columbia (Figure 7). Chaetomiaceae were highly abundant in all ITS libraries in the present study, comprising ~0.5, 3, and 9% of total in situ, ¹²-C and ¹³C-libraries. Both Sordariomycetes (Spearman's r = −0.37, p = 0.06) and Actinobacteria (r = −0.39, p = 0.05) were negatively correlated with C:N ratio, while Sordariomycetes were positively (r = 0.42, p = 0.03) and Actinobacteria negatively (r = −0.38, p = 0.05) correlated with total carbon in mineral layer soil (in situ libraries). The relative abundance of the aforementioned taxa did not significantly differ among the three sites (Blodgett, Brandy City, and Lowell Hill), while some taxa exhibited soil layer preferences (Table 2).

Indicator species analysis identified bacterial taxa with consistently increased relative abundance in soils from harvested plots across all pyrotag libraries (Figure 5; Supplementary Data 3). All of these taxa have members who are reported to be tolerant of heat, radiation and desiccation: Geodermatophilus (Montero-Calasanz et al., 2014; Sghaier et al., 2016), Sporichthya (Eppard et al., 1996; Babalola et al., 2009), Ramlibacter (De Luca et al., 2011), Flavisolibacter (Joo et al., 2015), Methylobacterium (Nogueira et al., 1998; Rokitko et al., 2003) and Segetibacter (Liu et al., 2014). The relative abundance of these groups (except Flavisolibacter and Sporichthyaceae) was also substantially increased in harvested plots from LTSP field sites in British Columbia (Figure S5).

A total of 234 bacterial ¹³C-enriched OTUs (enrOTUs) were identified, representing nine phyla. EnrOTUs classified as Actinomycetales, Armatimonadetes, Cytophagales, Myxococcales, Planctomycetes, Rhizobiales, Opitutaceae and Oxalobacteraceae were the most highly enriched (Table 2). Non-metric multidimensional scaling confirmed broad differences between ¹³C- and ¹²C-pyrotag libraries as well as distinct cellulolytic bacterial populations in each soil layer (Figure 3). EnrOTU from the organic-rich soil layer were mainly represented by previously known cellulose-degrading phyla, Cytophaga and Actinobacteria, while the mineral layer contained Betaproteobacteria and less characterized phyla such as Armatimonadetes (candidate division FBP and order FW68), Verrucomicrobia (classes Opitutae and Spartobacteria) and candidatus Saccharibacteria (formerly TM7). The delineation of fungal enrOTUs was less successful due to sparse overlap amongst OTUs in ITS libraries, which were typically dominated by few, highly abundant OTUs. This was illustrated by the poor separation of samples by NMDS (Figure S6) and the relatively small number of fungal enrOTUs identified (n = 16). These enrOTUs included unclassified Ascomycota and members of Agaricomycetes and Sordariomycetes (Table 2), while clustering in NMDS suggested the involvement of Dothideomycetes and a large proportion of unclassified ITS sequences. Ascomycota were major cellulose degraders under our experimental conditions as evidenced by the massive difference between ¹³C- vs. ¹²C- metagenomes in the proportion of reads classified as Ascomycota (an average of 10.6 and 0.8%, respectively; Figure S7).

Relative to the ¹²C-metagenome, unassembled ¹³C-metagenomes encoded double the number of glycosyl hydrolases (GH) and three-fold more GH families with reported endoglucanase activity. Five endoglucanase-containing families and lytic polysaccharide monooxygenases (AA9) were among the most enriched CAZy gene families in ¹³C-metagenomes (Figure 8). Lignin modifying enzymes, peroxidases (AA2) and iron reductase domains (AA8), were also highly enriched in ¹³C-metagenomes and were classified to the fungal order Sordariales, which contains the family Chaetomiaceae. The majority of differentially abundant GH genes were actinobacterial and fungal (Sordariales), while a lesser number were from Bacillales, Bacteroidales, Burkholderiales, Cytophagales, Opitutales, and Planctomycetes. Improved assemblies enabled the recovery of 10 taxonomically uniform MAGs of putatively cellulolytic bacteria (Figure 4; Table S3). The most complete were related to Myceliophthora thermophila (Ascomycota), Kitasatospora (Actinobacteria), Opitutaceae (Verrucomicrobia), Herbaspirillum (Betaproteobacteria), Chthoniobacter (Verrucomicrobia) and Caulobacteraceae (Alphaproteobacteria), though all likely represent, to an extent, a mixture of sub-populations.

Our results, based on multiple data types, demonstrate that timber harvesting can effect long-term changes in the composition of cellulolytic populations that reduce the rate of cellulose decomposition. This finding builds upon previous research that showed long-term impacts of harvesting on the composition and putative lignocellulolytic capacity of soil decomposers (Cardenas et al., 2015; Leung et al., 2016) by directly linking changes in composition with functional activity. Our conclusions differ from the large-scale meta-analysis of harvesting impacts on microbial activity by Holdena and Treseder (2013), since we did not observe a correlation between changes in activity and microbial biomass. This difference is not contradictory, rather, our results reveal the relative influence of community composition on activity when broad changes in the structure of decomposer populations occur. The most pronounced change in cellulose-degrading populations was an increase in the relative abundance of saprotrophic fungi, namely Chaetomiaceae, and a decrease in bacteria, notably Gram-positive Streptomycetaceae. These trends were robust across all datasets from both microcosm experiments and field samples and, though an effect on these specific groups has not been previously identified, similar trends in the relative abundance of saprotrophic fungi and Actinobacteria have been observed in soil from harvested forests after seven (Lewandowski et al., 2015), 15 (Hartmann et al., 2009, 2012) and fourty years (Chatterjee et al., 2008). In one case, populations of Sordariales, including the family Chaetomiaceae, were increased in logged forests in Southeast Asian tropical forests (McGuire et al., 2015). These communities were similarly impacted by forest fire (Xiang et al., 2014) and large-scale tree die back due to insect herbivory (Štursová et al., 2014), suggesting these trends may be broadly associated with forest disturbance.

Our observation that lower cellulolytic activity corresponded with increased relative abundance of fungi is at odds with the conventional view that fungi are the most effective decomposers of recalcitrant plant polymers in soil (Štursová et al., 2012). There are several cases in which soil properties, nutrient availability and litter quality influenced whether decomposition was predominantly fungal or bacterial (Jastrow et al., 2007; Güsewell and Gessner, 2009; Strickland and Rousk, 2010). Consistent with our results, Strickland et al. (2009) found that the relative abundance of Sordariomycetes was negatively correlated with net respiration during litter decomposition, while the reverse was true for Actinobacteria. A decrease in respiration activity could result from lower inherent activity of certain members of Sordariomycetes or, that their manner of decomposition affects the overall utilization of cellulose by other populations by either changing the quality of OM or producing inhibitory by-products. Whatever may be driving this ecological phenomenon, it is certainly of interest for understanding the effects of forest disturbance on nutrient cycling, as well as how community composition affects the carbon sequestration in soils.

The greatest impacts of harvesting on cellulolytic populations (in microcosms and in situ) were observed at the highest intensity of OM removal (OM3), while intermediate levels (OM1 and OM2) were generally indistinguishable from each other. These observations agree with the conclusions from previous studies characterizing the effects of OM removal on soil communities (Hartmann et al., 2012; Cardenas et al., 2015; Leung et al., 2016). The general lack of difference between OM1 and OM2, which differ by the retention of woody debris in OM1, suggests that the changes we observed most likely result from long-term exposure to stress-inducing environmental conditions, like heat and dryness, which were most pronounced in OM3. Five years following harvesting at the Californian sites, soil temperatures in harvested plots were between 5% (OM1) and 40% (OM3) higher than in unharvested plots, and soil moisture declined throughout the soil column in OM3 to a level deemed an erosion hazard (Paz, 2001). In general, near-surface soils, like those sampled in this study, experience significant abiotic changes in the interim between harvest and canopy closure, including higher average and maximum temperatures and lower moisture content, as well as greater fluctuations (Childs and Flint, 1987; Adams et al., 1991; Kranabetter and Chapman, 1999; Redding et al., 2003; Kulmala et al., 2014). The overall increase in stress-tolerant taxa in harvested soils was indicative of the influence of long-term changes in soil temperature and moisture regimes.

The expansion of cellulolytic members of Chaetomiaceae (dark-septate, thermophilic fungi), exemplified the general shift in populations adapted to harsher conditions. Chaetomiaceae are prevalent in hot, arid environments (Powell et al., 2012) and fire-prone forests (Rajulu et al., 2014) with several characterized thermophilic species, such as Myceliophthora thermophila (Berka et al., 2011). Similarly, the dramatic increase in harvested plots of non-cellulolytic taxa which possess remarkable stress-tolerance, such as Methylobacterium (Nogueira et al., 1998; Rokitko et al., 2003) and Geodermatophilus (Sghaier et al., 2016), is consistent with the influence of harsher environmental conditions. The decline in the abundance of cellulolytic Verrucomicrobia and Caulobacter in harvested plots may also be explained by drier soils in the decades post-harvesting. Verrucomicrobial populations have been positively correlated with soil moisture (Buckley and Schmidt, 2001), and Caulobacter species are generally known as aquatic organisms and have been found to respond rapidly to soil wetting (Fazi et al., 2008). The consistency of these trends across multiple datasets within our study and with previously published data from LTSP installations in British Columbia indicates the long-term selection pressures for stress-tolerant groups present post-harvesting.

For the most part, the retention of woody debris (OM1 relative to OM2) did not cause major differences in cellulolytic activity or composition of cellulolytic populations. An increase in the diversity of wood rot fungi can occur after harvesting when woody debris is retained (Brazee et al., 2014), but a similar trend was not supported by our data. One major point of difference between OM1 and OM2 was the relative abundance of fungi and bacteria, where OM1 and REF both had significantly more bacterial biomass than OM2 and OM3. There was some evidence to suggest differences in organic matter factored since both REF and OM1 had higher total carbon content as well as a slightly higher C:N ratio. Yet, these differences, for the most part, were not statistically significant and the correlations of major bacterial and fungal taxonomic groups (Actinobacteria and Sordariomycetes) to C:N (−/−, respectively) and total carbon (±) were at odds with the trends in fungal and bacterial biomass among treatments. While our study was not designed to discern between the effects of OM retention on environmental conditions vs. on OM quality and quantity, the clearest impact of OM retention in our results was to mitigate changes in the relative abundance of stress-tolerant taxa.

Changes in cellulolytic populations differed in several ways from previous observations of the long-term effects of harvesting on whole soil communities. Firstly, a decrease in the ratio of Basidiomycota to Ascomycota is a common characteristic of soil fungal communities in the years and decades following timber harvesting (Bader et al., 1995; Hartmann et al., 2012; Štursová et al., 2014; McGuire et al., 2015) and following forest fire (Holdenb et al., 2013; Buscardo et al., 2015). Accordingly, we found a significant decrease in the ratio in harvested plots from field samples (i.e., in situ), but not in cellulolytic populations (i.e., ¹³C-pyrotag libraries). The lack of difference among cellulolytic populations may simply reflect the predominance of Ascomycota as cellulose degraders, which overshadowed any negative impacts of harvesting on cellulolytic Basidiomycota. In general, the predominance of cellulolytic Ascomycota in forest soils supports the hypothesis that the generally observed shift in the ratio of Basidiomycota and Ascomycota reflects a post-harvesting shift to a saprotroph-dominated system (Hartmann et al., 2012). Notably, the relative abundance of one of the two cellulolytic groups of Basidiomycota identified, Clitopilus spp., increased in OM3, suggesting some Basidiomycota may thrive in post-harvest conditions. Secondly, the diversity of cellulolytic populations was not impacted by harvesting, contrary to a similar study of the effects of prescribed burning which found a significant decrease in the diversity of cellulolytic taxa (Bastias et al., 2009). In general, soil microbial diversity does not differ between primary to secondary forests (Lauber et al., 2008; Paula et al., 2014; McGuire et al., 2015; Oliver et al., 2015), but is commonly reduced by the conversion of forest to agricultural land (Rodrigues et al., 2013). Thus, our findings for cellulolytic populations were in broad agreement with previous characterizations, suggesting diversity is less indicative of harvesting effects compared to shifts in community structure, in our case, toward stress-tolerant taxa.

The majority of taxa identified by SIP (~75%) had previously documented cellulolytic activity, including well-characterized groups, such as Actinobacteria, Bacteroidetes, Cytophaga, Myxococcales, and Sordariomycetes. Though the potential of non-cellulolytic taxa incorporating sufficient ¹³C-label via cross-feeding to be designated as cellulolytic cannot be ruled out in SIP experiments, this level of agreement with culture-dependent characterizations supports the effectiveness of SIP as a culture-independent method to link sequence data with function. The remaining taxa, not previously known to degrade cellulose, were largely associated with mineral layer soil, an atypical sample source for studies of cellulose-degradation. For example, well-characterized cellulolytic taxa, like Cytophaga and Actinobacteria, were associated with organic layer soils, while mineral layer-associated cellulolytic taxa belonged to candidate division FBP (Armatimonadetes), with no representative genome or isolate, and members of the ubiquitous, yet poorly characterized phylum, candidatus Saccharibacter (formerly TM7), of which we recovered a partial MAG (~0.4 Mb). Phyla with relatively few cultured representatives, such as Armatimonadetes, Verrucomicrobia, and Planctomycetes, were also mineral layer-associated, each possessing at least one representative characterized to degrade cellulose (Sangwan et al., 2004; Dedysh and Kulichevskaya, 2013; Lee et al., 2014). In terms of cellulolytic activity, microbial biomass in mineral soils exhibited higher activity per unit total biomass than organic layer soil, and, in mineral layer soil, ¹³C was primarily assimilated by bacteria. Thus, this study revealed mineral soil to possess active cellulolytic populations, distinct from those in the organic layer soils, comprised of poorly characterized cellulolytic taxa.

SIP-based designations of cellulolytic fungi also matched previously characterized cellulolytic taxa, such as Humicola, Clitopilus, and members of Chaetomiaceae. The putative assignment of a member of Sebacinaceae was novel, given members of this family are better known for their ectomycorrhizal associations. Yet, a few saprobic species of Sebacinaceae have been described, though they have not previously been reported to be cellulolytic (Oberwinkler et al., 2013; Weiß et al., 2016). The 46-Mb fungal MAG was the first eukaryotic genome of its size recovered from shotgun metagenomic data and, a size that resembles that of a Myceliophthora thermophila draft genome (~38.7 Mb; Berka et al., 2011). Our MAG encodes a number of CAZy families with endoglucanases, lytic polysaccharide monooxygenases (AA9), iron reductase domains (AA8) and peroxidases (AA2), which suggest a role in both cellulose and lignin decomposition that mirrors cultured representatives of M. thermophila, which can completely degrade lignocellulose and encode a diverse array of thermostable CAZymes (Babot et al., 2011; Berka et al., 2011). The recovery of multiple MAGs, including such a large one, demonstrates the power of combining SIP and shotgun metagenomics to investigate the function of taxa in highly diverse soil communities where metagenome assembly is otherwise typically poor.

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.