We processed the forest floor in a blender to pass a 2 mm mesh screen and then pulverized a subsample of the blended material in a ball mill (SPEX Sample Prep LLC, Metuchen, NJ, United States). We oven-dried the pulverized forest floor material at 65°C prior to determination of total C and N. We used the sieved mineral soil for determination of total C, N, pyrogenic-C (PyC), NO₃-N, NH₄-N, extractable P, and pH. We oven-dried the mineral soil sample to be used for C and N analysis at 105°C and then pulverized the subsamples as described above. We analyzed one forest floor sample and one mineral soil sample from each subplot for total C and N concentrations on a dry combustion elemental analyzer (Costech Analytical Technologies Inc., Valencia, CA, United States), using atropine as the quantification standard. We quantified mineral soil PyC concentrations by digesting 0.5 g mineral soil from each subplot in 10 mL 1.0 M HNO₃ + 20 mL 30% H₂O₂ at 100°C for 16 h (Kurth et al., 2006) and analyzing post-digested soils for residual C and N.

We extracted a 10 g sample of fresh mineral soil from each subplot for NH₄⁺-N and NO₃^–-N in 50 mL 2.0 M KCl for 1 h. We filtered the resulting soil extracts through 2.5 μm pore-size filter paper (GE Healthcare UK Limited, Little Chalfont, Buckinghamshire, United Kingdom). We determined the extract NH₄⁺ concentrations spectrophotometrically by reacting with ammonia salicylate and ammonia cyanurate (Sinsabaugh et al., 2000), and measuring absorbance at 595 nm (BioTek Elx800, BioTek Instruments Inc., Winooski, VT, United States). We determined the concentrations of NO₃^– in the extracts spectrophotometrically by reacting the extracts with vanadium (III) chloride, sulfanilamide, and N-(1-naphthyl)-ethylenediamine dihydrochloride and measuring absorbance at 540 nm (Doane and Horwáth, 2003). We extracted 2.5 g subsamples of fresh mineral soils for phosphorus (P) in 40 mL 0.5 M NaHCO₃ and determined P concentration using the Olsen method (Olsen et al., 1954). We measured mineral soil pH of a 1:2 (w:v) soil slurry with a benchtop pH meter (Oakton pH 700, Oakton Instruments, Vernon Hills, IL, United States).

We isolated the light particulate organic matter (LPOM) fraction from mineral soil to determine SOM composition via diffuse reflectance Fourier transform infrared (DRIFT) spectroscopy. To collect LPOM, we first dispersed mineral soil in DI water via sonication. We then passed the soil solution through a 53 μm sieve, followed by isolation of LPOM via flotation in 1.7 g cm^–3 sodium iodide solution. We then pulverized the LPOM fraction using mortar and pestle. We scanned undiluted LPOM samples in aluminum 96 well-plates in diffuse reflectance mode using a mid-infrared spectrometer (Vertex 70, Bruker Scientific LLC, Billerica, MD, United States). We scanned from 4,000 to 400 cm^–1 with a resolution of 4 cm^–1 and 64 scans per spectrum. Background spectra from empty wells were collected with every plate. We background-subtracted and baseline corrected the resulting spectra in OPUS 7.5 (Bruker Scientific) and calculated the relative areas of peaks indicative of aliphatic (2,930–2,870 cm^–1), carbohydrate (1,160–1,040 cm^–1), aromatic (1,590–1,570 cm^–1, 1,550–1,500 cm^–1, 975–700 cm^–1), and amide (1,450–1,400 cm^–1, 1,320–1,220 cm^–1) compounds (Calderón et al., 2011). We considered the sum of aliphatic and carbohydrate peak areas as a proxy for labile C.

We incubated soils to determine the size and kinetic rates of the C_a and C_s pools and potential soil C flux rates (Paul et al., 2001). We adjusted a 30 g sample of fresh mineral soil from each subplot to 40% water filled pore space (WFPS) in 120 mL specimen cups. The specimen cups were placed in 1 L glass jars, and the soils were incubated in the dark for 300 days at ambient temperature (∼23°C) with biweekly adjustments of soil moisture to 40% WFPS. We measured CO₂ evolution on days 10, 14, 28, 42, 58, 90, and every 30 days thereafter until day 300. Prior to each measurement event, we flushed the jars to ambient CO₂ concentrations, then tightly sealed the jars for 24–48 h before sampling a 1 mL gas aliquot through septa fitted to the jar lids. We measured CO₂ concentration of the aliquot using an infrared gas analyzer (LI-COR Inc., Lincoln, NE, United States), and calculated CO₂-C efflux as the difference in CO₂-C concentrations between the soil containing jars and blank jars that contained only a specimen cup and DI water.

We extracted DNA from 0.25 g of mineral soil from each of the 40 subplots using the MoBIO PowerSoil DNA isolation kit (MoBIO laboratories, Carlsbad, CA, United States), according to the manufacturer’s instructions. The extracted DNA was amplified using PCR with barcoded 515f/806r universal primers of the V4 region of the 16S-rRNA gene (Caporaso et al., 2010). Samples were sequenced on the Illumina MiSeq platform using the v2 500 cycle Reagent Kit at the Michigan State University Genomics Core. We processed DNA sequences using the QIIME2 bioinformatics pipeline (Bolyen et al., 2019). We denoised, merged forward and reverse reads, and removed chimeras using the q2-DADA2 plugin (Callahan et al., 2016). After sequence processing there were 20,815–67,314 sequences per sample (mean = 36,766). We did not rarefy our samples because rarefying decreases statistical power and is not recommended when library sizes vary by less than 10-fold (McMurdie and Holmes, 2014; Weiss et al., 2017). We inferred phylogenetic trees by applying MAFFT multiple sequence alignment (Katoh and Standley, 2013) and FastTree 2 (Price et al., 2010) using the q2-phylogeny plugin. With the q2-diversity plugin, we used the phylogenetic tree to calculate OTU richness, Faith’s phylogenetic diversity (Faith, 1992), and a weighted UniFrac distance matrix (Lozupone and Knight, 2005). We used the q2-feature-classifier plugin (Bokulich et al., 2018) with the classify_sklearn action (Pedregosa et al., 2011) to classify taxonomic composition of our samples employing a Naïve Bayes classifier trained on the SILVA SSURef database version 132 using a 99% similarity threshold (Quast et al., 2013). We filtered OTUs associated with Archaea, Eukaryotes, mitochondria, and chloroplasts. We calculated oligotroph-to-copiotroph ratio at the phylum level as the ratio of the sum of relative abundances of all taxa classified within the phyla Acidobacteria and Verrucomicrobia to the sum of the relative abundances of all taxa classified within the phyla Actinobacteria, Firmicutes, and Bacteroidetes (Fierer et al., 2007, 2012; Ramirez et al., 2012). We predicted 16S rRNA gene copy numbers using the ribosomal RNA operon database (rrnDB) (Stoddard et al., 2015). For each OTU, we assigned 16S gene copy number as the genus mean; if a genus was not represented in rrnDB we assigned the family mean and so on up taxonomic levels until all OTUs were assigned a 16S gene copy number value. We then calculated the abundance-weighted mean 16S gene copy number for each bacterial community as another indicator of community-scale copiotrophy (Nemergut et al., 2016). This approach of predicting 16S gene copy numbers based on databases has frequently been used as an indicator of community-scale life-strategy (Nemergut et al., 2016; Roller et al., 2016; Kurm et al., 2017; Wu et al., 2017; Li et al., 2021), but these estimates can be biased based on reference genomes in the database and should thus be interpreted with caution (Louca et al., 2018).

We predicted metagenomic functional potential of soil bacterial communities using the PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) software (Douglas et al., 2019). PICRUSt2 uses HMMER software (hmmer.org) to perform multiple sequence alignment via hidden Markov models, places reads in a reference tree using EPA-ng (Barbera et al., 2018), and outputs a tree file with Gappa (Czech et al., 2020). PICRUSt2 then performs hidden-state prediction of gene families (Louca and Doebeli, 2017), and uses MinPath (Ye and Doak, 2009) to infer the relative abundance of MetaCyc metabolic pathways (Caspi et al., 2018). MetaCyc is a metabolic pathway database in which pathways are hierarchically grouped according to metabolic function. We focused our analysis on degradation pathway groups related to C cycling, including carbohydrate degradation, alcohol degradation, amine degradation, amino acid degradation, aromatic compound degradation, carboxylate degradation, fatty acid and lipid degradation, nucleoside and nucleotide degradation, and secondary metabolite degradation pathways. PICRUSt2 calculates weighted nearest sequenced taxon index (NSTI), which is an indicator of the average phylogenetic distance of bacteria within an OTU table compared to genomes in the reference database, and is correlated with accuracy of metagenome predictions (Langille et al., 2013). The mean weighted NSTI for our samples was 0.12 (SD = 0.04), which is similar or better to that in other soil studies (Canarini et al., 2021; Toole et al., 2021; Wang et al., 2021), and indicates that metagenomic functional profiles can be accurately predicted in our dataset.

We performed all statistical analysis in the R statistical computing environment (R Core Team, 2019). We used linear mixed models (R Package: nlme; Pinheiro et al., 2019) to determine the impacts of wildfire occurrence and burn severity (dNBR) on live and dead tree basal area, total C and N, C:N ratio, PyC, NH₄⁺-N, NO₃^–-N, TIN, extractable P, pH, gravimetric soil moisture, forest floor mass, and relative DRIFT peak areas. To account for potential legacy effects of soil properties related to topography, all models initially included elevation, slope, and aspect as covariates. Topographic variables that were not significant at α < 0.05 were sequentially removed from the models. All models included plot-identifier as a random effect, and the NH₄⁺-N, NO₃^–-N, and TIN models included total N as a covariate. We determined that inorganic N variables did not meet assumptions of normality, so we normalized these variables using Box-Cox transformations (Box and Cox, 1964). We determined whether wildfire occurrence or dNBR was a better predictor of differences by fitting separate models that included either fire occurrence (burned or unburned) or burn severity as an explanatory variable and comparing model AIC values. We selected fire occurrence or severity as the better predictor if the model containing that variable was ≥ 2.0 AIC points lower than the alternative model.

We assessed the impacts of fire and severity on soil C flux rates and cumulative soil C flux over our soil incubation using the following fixed effects portion of linear mixed models:

where β_i represents the model-fitted intercept or slope coefficient and fire variable is either fire occurrence or burn severity. Each model also initially included all topographic variables and a plot-identifier random effect. Non-significant topographic variables were sequentially removed, and selection of fire occurrence or severity as a more suitable predictor was performed using AIC values as described above.

We assessed the size and mineralization rates of the C_a and C_s pools using non-linear mixed models. We assessed the sizes of the C_a and C_s pools and their mineralization rate constants using the following fixed effects portion a two-compartment first-order kinetics model (Kuzyakov, 2011):

where ^RateCO₂ is the soil CO₂ efflux rate at each measurement event, C_a is the size of the active C pool, k_a is its mineralization coefficient, C_s is the size of the non-active C pool and k_s is its mineralization rate coefficient. In this model, C_s is constrained to be C_SOC – C_rma, where C_SOC is total soil organic C content.

We determined whether fire occurrence or severity affected soil C pools by first fitting a global non-linear model that included all 40 incubation replicates. We performed model selection procedures for random effects in non-linear models as described by Pinheiro and Bates (2000), which resulted in the inclusion of a plot identifier as a random effect associated with C_a and k_s, and no random effect associated with k_a. We assessed how wildfire impacted the size or kinetic rates of the C_a and C_s pools by adding fire occurrence or severity as a covariate to the non-linear models. Selection of fire occurrence or severity as a more suitable predictor variable was performed using AIC values as previously described. In addition to the global model, we fit separate non-linear models for each incubation replicate to obtain independent estimates of C_a, k_a, and k_s for use in elastic net regularization analysis and structural equation modeling. Five incubation replicates (two replicates from separate unburned plots, and one each from plots with dNBR values of 108, 350, and 478) did not converge during non-linear regression due to consistently high respiration rates throughout the incubation and were excluded from further analysis. After exclusion, each of the ten plots still had C pool parameter estimates from at least three subplots.

We assessed the impacts of wildfire and burn severity on Faith’s phylogenetic diversity, OTU richness, Shannon’s diversity, Pielou’s evenness, O:C, weighed 16S rRNA gene copy number, and the relative abundance of the ten most abundant bacterial phyla using linear mixed models and AIC model-selection as described above. The ten most abundant bacterial phyla represented all phyla that exhibited > 0.5% relative abundance and 94.1% of the total abundance. We used multivariate statistical approaches (R Package: vegan; Oksanen et al., 2019) to determine the impacts of fire occurrence, burn severity, and soil properties on bacterial community structure. Using a weighted UniFrac distance matrix of the bacterial community, we employed principle coordinates analysis (PCoA) to determine the relationships between plant and soil characteristics and bacterial communities. Additionally, we performed PCoA on a Bray-Curtis distance matrix of imputed MetaCyc metabolic pathways, hierarchically grouped into common C-degradation pathways.

We identified positive and negative fire responsive bacterial OTUs using species indicator analysis with fire occurrence as the explanatory variable, relative abundance of OTUs as response variables, and point biserial correlation as the output statistic (R Package: indicspecies; De Caceres and Legendre, 2009). We assessed which positive fire-responder OTUs were also positively related to severity by performing linear mixed-model analysis on the relative abundance of each fire-responsive OTU with dNBR as the explanatory variable and unburned plots excluded from the analysis.

We determined properties linked to the size and kinetic rates of soil C pools using structural equation modeling. We constructed an initial SEM meta-model composed of linear mixed models that included live tree basal area, dead tree basal area, forest floor mass, total C, PyC, total N, TIN, extractable P, pH, and soil moisture as response variables (R Package: piecewiseSEM; Lefcheck, 2016). We assessed relationships of soil properties, live and dead tree basal area, and topography with C pools using separate SEMs for each C pool parameter (i.e., C_a, k_a, k_s). Additional methodological details and results of our SEMs are provided in Supplementary Note 1.

We used correlation analyses to assess relationships between log₂ bacterial abundance at the levels of phylum and genus with C pools, soil TIN and P, pH, and relative DRIFT peak areas associated aromatic and labile compounds in the LPOM fraction. We selected these soil variables because they indicate soil nutrient availability and the amount of labile vs. recalcitrant C, and are therefore associated with bacterial copiotrophic versus oligotrophic life-strategies. We performed univariate correlations between all bacterial phyla and the selected soil variables. At the genus level, we performed correlations within phyla that our linear mixed-models indicated were related to burn severity and that previous research indicates exhibit either copiotrophic or oligotrophic life strategy (i.e., Acidobacteria, Actinobacteria, Bacteroidetes, Proteobacteria, and Verrucomicrobia). All correlations were performed using Spearman’s rank correlation to account for non-normality of bacterial abundance data and we considered correlations significant at α = 0.05.

We performed elastic net regularization (R Package: glmnet; Friedman et al., 2010) to identify bacteria phyla that are predictive of C pool parameters. Elastic net regularization selects the best set of predictors for a given response by using LASSO and ridge regression, but does not assign significance values to selected predictors; elastic net regularization is frequently used with sparse datasets when the number of potential predictors is larger than the number of observations, as is often the case for genomic data (Friedman et al., 2010; Wagg et al., 2019). We used linear mixed-models to determine relationships between C degradation metabolic pathways and C pool parameters.

Our AIC-based model selection indicated that fire occurrence had more explanatory power for models assessing total soil C, NH₄⁺-N and TIN, whereas severity had more explanatory power for dead tree basal area, forest floor mass, total N, pH, and soil moisture (Table 1). Neither explanatory variable had more power for live tree basal area, C:N, NO₃^–-N, extractable P, or relative DRIFT peak areas.

Concentrations of NH₄⁺-N, NO₃^–-N, and TIN were higher in burned stands than unburned stands (Table 1). Additionally, total N was a significant covariate for NH₄⁺-N, NO₃^–-N, and TIN, all of which were all positively related to total N (p < 0.001, p = 0.041, and p < 0.001, respectively). Dead tree basal area increased with severity (p = 0.005), whereas live tree basal area was not related to severity. Forest floor mass was negatively associated with severity (p = 0.002). Mineral soil gravimetric moisture was negatively related to severity (p = 0.004). Soil pH increased in association with severity (p = 0.007). Mineral soil total C and N, C:N, PyC, extractable P, and relative DRIFT peak areas were not significantly affected by fire occurrence or burn severity. Our SEM meta-model revealed relationships between severity and soil properties that were not captured by linear mixed-models (Supplementary Note 1 and Supplementary Figure 2).

Fire occurrence and severity negatively impacted Faith’s phylogenetic diversity, with AIC values indicating no preference for fire occurrence or severity-based models (Figures 2A,B). AIC values indicated that the severity model better explained OTU richness, which was negatively correlated with severity (Figures 2C,D). Fire occurrence had more explanatory power for changes in phyla-level oligotroph-to-copiotroph ratios, which decreased from 1.08 in unburned soils to 0.48 in burned soils (Figures 2E,F). Fire occurrence and severity were positively associated with 16S rRNA gene copy number, with AIC values indicating no preference for either explanatory variable (Figures 2G,H).

We performed PCoA on a weighted UniFrac distance matrix and found that the first axis explained 54.9% of variation, and the second axis explained an additional 8.7% (Figure 3). Fire occurrence (r² = 0.24, p < 0.001) and severity (r² = 0.57, p < 0.001) were both significantly correlated with the PCoA ordination. Several soil properties were also significantly correlated with the ordination, including NH₄-N concentration (p = 0.002), TIN concentration (p = 0.008), P concentration (p = 0.009), pH (p < 0.001), soil moisture (p = 0.010), forest floor mass (p = 0.003), and labile C peak areas (p = 0.020) and aromatic C peak areas (p = 0.024) in the LPOM fraction.

We assessed the relationship between wildfire occurrence and severity and the relative abundances of the ten most abundant bacterial groups using univariate linear mixed models. AIC values indicated that the fire occurrence models had more explanatory power for assessing differences in Bacteroidetes, Acidobacteria, Verrucomicrobia, and Planctomycetes relative abundance. The relative abundances of Bacteroidetes, Actinobacteria, and Firmicutes were higher in burned areas than unburned areas (Figure 4). The relative abundances of Acidobacteria, Verrucomicrobia, and Planctomycetes relative abundance were lower in burned areas than unburned areas. There was no significant effect of fire occurrence or burn severity on α-Proteobacteria,γ-Proteobacteria, δ-proteobacteria or Gemmatimonadetes relative abundance.

Indicator species analysis identified 53 OTUs as indicators of burned soils (i.e., positive fire responders) and 74 OTUs as indicators of unburned soils (i.e., negative fire responders) (Supplementary Tables 1, 2). The positive fire responders most commonly belonged to the Actinobacteria and Bacteroidetes phyla, as well as to the α-proteobacteria class, which, respectively, accounted for 28.6, 23.2, and 26.8% of the positive responder OTUs. The seven OTUs that exhibited the strongest positive response to fire (point biserial correlation > 0.60) came from the genera Massilia (γ-proteobacteria), Roseomonas (α-proteobacteria), Segetibacter (Bacteroidetes; 2 OTUs), Blastococcus (Actinobacteria), unclassified Micrococcaceae genus (Actinobacteria), and unclassified Burkholderiaceae genus (γ-proteobacteria). The negative fire responders most frequently belonged to the Planctomycetes phylum, which accounted for 14.7% of these responders, and to the α-proteobacteria and γ-proteobacteria classes, which accounted for 18.7 and 17.3%, respectively. The five OTUs that exhibited the strongest negative response came from the genera Cytophaga (Bacteroidetes), IS-44 (γ-proteobacteria), Mycobacterium (Actinobacteria), uncultured Elsteraceae genus (α-proteobacteria), and uncultured Gemmataceae genus (Planctomycetes). OTUs identified as positive fire responders had a significantly higher mean 16S gene copy number than negative fire responders (3.41 vs. 2.64; p = 0.020).

We identified 15 OTUs as positively associated with severity based on linear mixed models (Supplementary Figure 4). Of the 15 OTUs, six were from the Bacteroidetes phylum, two were from Actinobacteria, and one was from Verrucomicrobia. Four severity responders were from the α-proteobacteria class, one was from the γ-proteobacteria class, and one was from the δ-proteobacteria class. The abundances of all of the severity-associated OTUs were positively correlated with either TIN or P (data not shown).

Using PICRUSt2, we identified 135 MetaCyc metabolic pathways associated with C-degradation functions. PCoA on a Bray-Curtis matrix of the MetaCyc pathways grouped into common C-degradation functions indicated that the first axis explained 45.1% of variation and the second axis explained 22.1% (Figure 5). Fire occurrence (r² = 0.16, p < 0.001) and severity (r² = 0.41, p < 0.001) were significantly correlated with the ordination. Additionally, several soil properties were significantly correlated with the ordination, including total C (p = 0.009), C:N ratio (p = 0.044), NH₄-N (p = 0.008), TIN (p = 0.047), extractable P (p = 0.002), pH (p = 0.010), and forest floor mass (p = 0.007).

AIC values indicated that severity had more explanatory power for differences in carbohydrate degradation, alcohol degradation, amine and polyamine degradation, carboxylate degradation, and nucleotide and nucleoside degradation. Severity positively impacted the relative abundance of imputed pathways associated with carbohydrate degradation, alcohol degradation, amine and polyamine degradation, carboxylate degradation, and secondary metabolite degradation (Figure 6). The relative abundance of amino acid degradation pathways was higher in burned areas than unburned areas (p = 0.021). The relative abundances of aromatic compound degradation pathways and fatty acid and lipid degradation pathways were not related to fire occurrence or severity.