Wildfires are increasing in frequency and severity with changes in the climate. Increases in fire activity may lower soil carbon (C) stocks if photosynthetic C inputs are outpaced by microbial decomposition in burned ecosystems. However, fires may preferentially deplete particulate organic carbon (POC, which is more susceptible to microbial decomposition) over mineral-associated organic carbon (MAOC, which is protected from decomposers), potentially slowing microbial C losses after wildfires. Yet it remains unclear how plants, microorganisms, and soil organic matter pools interact to control the fate of C after wildfires. To assess how wildfires affect the persistence of soil C, we measured POC, MAOC, pyrogenic organic C, plant cover, extracellular enzyme activity (EEA), soil microbial abundance, and microbial community composition 17 days, 1, 3, and 4 years after the Holy Fire burned 94 km2 of a fire-adapted chaparral. We found that the fire immediately decreased POC by 50% (from 50.9 ± 20.5 to 25.5 ± 5.67 µg C g-1) and MAOC by 33% (from 9.25 ± 0.87 to 6.31 ± 0.92 µg C g-1), consistent with our hypothesis that MAOC would be less vulnerable to loss than POC. POC further decreased by another 38% by 1 year post-fire, likely as a result of increased microbial access to this C given increases in microbial abundance and EEA. Between one and four years after the fire, the percent cover of the dominant shrub (Arctostaphylos glandulosa) increased from 3.9 ± 1.6 % to 16 ± 5.4 % (compared to 58 ± 4.6 % in unburned control plots), marking the end of the decreasing trend in soil C. However, soil C stocks did not increase during this time, likely because plant C inputs did not meet microbial demand for C, a finding consistent with the enrichment in soil δ13C as microbes would have preferentially metabolized δ12C. As global change increases wildfire frequency and severity, C volatilization and decomposition may outpace the recovery of plant C inputs during the first four years after fire, slowing the recovery of soil C.

Site description

This study was conducted within the Holy Fire burn scar in the Cleveland National Forest. The Holy Fire burned 94 km2 of manzanita-dominated chaparral shrubland between August 6 and September 13, 2018. Seventeen days after the fire was extinguished, we established 9 plots; 6 of these plots were within the burn scar of the Holy Fire, while the other 3 served as control unburned plots (Figure S1; Pulido-Chavez et al. 2022). Each plot (~25 m2) consisted of four 1 m2 subplots located 5 m from the center of the plot in each cardinal direction (i.e., N, S, E, W; Figure S1). All nine plots had similar aspects, slope, elevation, and pre-fire vegetation dominated by manzanita (Arctostaphylos glandulosa) and chamise (Adenostoma fasciculatum) shrubs.

The climate at the site is Mediterranean with hot, dry summers (average 9.2 mm precipitation per month since 1990) and cool, wet winters (average 101 mm precipitation per month since 1990). Annual temperature averages 17 °C, and total annual precipitation averages 668 mm. Soils in all the control plots are classified as Lithic Haploxerolls and mapped in the Friant series, as are soils in half of the plots within the burn scar. Because of steep slopes, fire burn patterns, and limited accessibility to the burn scar, the remaining burned plots were established on Typic Xerorthents, mapped in the Cieneba series. Soils from both series are sandy and gravelly loams.

Soil Sampling

 At the first time point, we measured and averaged ash depth at 3 points within the 1 m2 subplots as an index of soil burn severity (Pulido-Chavez et al. 2022). We sampled soils from all 4 subplots within each of the 9 plots starting 17 days after the fire was extinguished, and returning 371, 1092, and 1460 days later (roughly one, three, and four years post-fire. All samples were collected in the fall to minimize the effect of seasonality on SOC. Soil samples from the A horizon (0–10 cm depth) were collected using an ethanol-sterilized bulb planter after removing the litter (unburned plots) and ash layers (burned plots) from the surface. Soils were transported back to the lab within hours of collection, where they were stored overnight at 4 °C. Soils were sieved to 2 mm within 24 hours of sampling. One subsample (~5-20 g) was frozen at -80 °C for molecular analyses, with the remaining soil left to air dry for physiochemical analyses. Soil pH was measured in a slurry of 2:1 water to dry soil. We only measured soil C concentrations and extracellular enzyme activity from soils collected from the N and S subplots. While microbial communities were measured in all four subplots (Pulido-Chavez et al. 2022), we only present microbial data from samples that were also analyzed for soil C.

Soil carbon pools

We used density fractionation to separate soil organic matter into particulate organic matter (POM), mineral-associated organic matter (MAOM), and sand-associated organic matter from the north and south subplots of each plot at each time point. Because little C was detected in the sand fraction (mean C in sand = 2.5 µg-C g-1), and MAOM plays a stronger role in stabilizing C (Poeplau et al., 2018; Püspök et al., 2022), we do not present sand-associated organic matter here. To fractionate organic matter, 5 g oven-dried soil (60 °C) was mixed with 1.65 g cm-3 sodium polytungstate (SPT) and dispersed using 200 J ml-1 ultrasonic energy at 60 W. The soil solution was then centrifuged at 2500 g for 60 minutes. The floating fraction was aspirated onto a 0.45 µm glass fiber filter for vacuum filtration and rinsing with deionized (DI) water. The remaining pellet was then washed with DI water and passed through a sieve to separate sand (> 53 µm) from the silt and clay (< 53 µm) fractions. The three soil fractions from each sample were then oven-dried (60 °C) for 24 hours, weighed, and finely ground before being analyzed for soil C and N content on an elemental analyzer (Flash EA 1112 CN; Thermo Fisher Scientific INC.). On average, we recovered 98 ± 4.9 % of soil mass, 117 ± 37 % of soil N, and 95 ± 25 % of soil C. We also measured δ13C of POC and MAOC to assess if plant C inputs and microbial respiration affected soil δ13C. We measured δ13C values from samples collected at the first (17 days, to determine a fire effect) and last (four years, to determine if values returned to pre-fire) timepoints at the UC Riverside Facility for Isotope Ratio Mass Spectrometry (FIRMS; https://ccb.ucr.edu/facilities/firms). 

To estimate the concentration of pyrolyzed organic matter (PyOM) in soils, we digested soil samples with nitric acid and hydrogen peroxide (Kurth et al., 2006; Licata & Sanford, 2012; Pingree et al., 2012). Soil samples (1 g) were oven dried (60 °C) before being mixed with 10 mL of 1 M nitric acid and 20 mL of 30% hydrogen peroxide in a 100 mL digestion tube. The solution was mixed by swirling the tube and allowed to sit at room temperature for at least 5 minutes. The digestion tube was then incubated at 100 °C for 16 hours, after which the samples were cooled before filtering (Whatman #2). The soil that remained on the filter paper was then dried before being analyzed for soil C content. Carbon in PyOM (PyOM-C) was calculated as the percent C times the mass of the digested soil sample divided by the weight of the initial undigested soil.

Soil potential extracellular enzyme activity and decomposition rates

We measured the potential activity of soil extracellular hydrolytic enzymes using substrates tagged with fluorescing methylum-belliferone (MUB) or 7-amido-4-methyl-coumarin (MC) (Marx et al., 2001). Activity of the starch-degrading enzyme α-glucosidase (AG) was measured using 4-MUB-α-d-glucopyranoside. Activity of the cellulose-degrading enzyme β-glucosidase (BG) was measured using 4-MUB ß-d-glucopyranoside. Activity of the cellulose-degrading enzyme 1,4-β-cellobiohydrolase (CBH) was measured using 4 MUB-β-D-cellobioside. Activity of the chitin-degrading and N-acquisition enzyme N-acetyl-ß-D-glucosaminidase (NAG) was measured using 4-MUB N-acetyl-ß-D-glucosaminide. Activity of the peptide-degrading and N-acquisition enzyme l-Leucine aminopeptidase (LAP) was measured using L-Leucine-7-amido-4-methylcoumarin. Finally, the activity of the organic P-degrading and P-acquisition enzyme Phosphomonoesterase (PHO) was measured using 4-MUB-phosphate.

 All enzyme assays were carried out in 96-well plates consisting of 100 µL substrate and 100 µL soil slurry. Soil slurries were made by blending soil in a universal buffer solution adjusted to pH 6.5 (1:125 soil weight to volume). For each substrate, we measured the background fluorescence of soils and substrate, and the quenching of MUB by soils, and used standard curves of MUB or MC to calculate the rate of substrate hydrolyzed. Fluorescence was measured at 360 nm excitation and 460 nm emission in a plate reader (Tecan, Männedorf, Switzerland) five times over a 3-h period following methods adapted from Saiya-Cork et al (2002) and McLaren et al (2017). We also measured the activity of the lignin-degrading enzymes phenol oxidase (Phenol) and peroxidase (Perox) using L-3,4-dihydroxyphenylalanine substrate. Color absorbance was measured at 460 nm in a plate reader after a 24-h incubation.

We measured leaf litter decomposition rates using a litter mesh bag approach in the second year post-fire. First, we collected recently senesced leaf litter from Arctostaphylos glandulosa shrubs adjacent to our control plots. This litter was brought to the lab, where it was shredded using a blade grinder to represent degradation via macro consumers. The ground litter (3 g) was then sealed within a pouch made from nylon screening (0.5 mm mesh size) and gamma irradiated at 12,000-14,000 Gy for 2.5 days at the University of California-Irvine, Dept of Radiation Oncology, to reduce native bacteria and fungi as previously established (Glassman et al. 2018). Approximately one month after collecting the litter, sterile decomposition bags were placed inside a cage made of wire (1 cm mesh size) to protect the mesh bags from destruction by insects and larger animals. A total of 3 litter bags per plot were randomly fixed to the ground using wire stakes on December 17, 2019. The first litter bag was collected 5 months after deployment, the second was collected at 13 months, and the third at 2 years. In the lab, the leaf litter was removed from the litter bags, dried overnight at 65°C, and reweighed to calculate mass loss. Decomposition rates were calculated as the percent of litter mass remaining after the field incubation.

2.5 Shrub and Grass Cover

We measured the cover of all plants in each of the 1m2 subplots at 286, 376, 612, 718, 966, and 1314 days after the fire was extinguished. Plant cover was measured by visually estimating vegetative cover for all species contained within each 1 m2 subplot. Plant cover was allowed to go over 100% for a given subplot due to overlapping plant canopies. We also measured the δ13C of leaves collected from dominant shrubs in burned and unburned plots. Within burned plots, leaves were collected from individual shrubs that had resprouted vegetatively or grown from seed. We attempted to sample six shrubs of each species from each type (adult, resprouter, or seedling), but only found five resprouted Adenostoma fasciculatum shrubs and three seeded Adenostoma fasciculatum shrubs. 

DNA extraction, sequencing, qPCR, and bioinformatics

We extracted DNA from frozen (-80° C) soils to measure the abundance, alpha diversity, and beta diversity of bacteria and fungi. For the 17-day, one-year, and three-year samples, we extracted DNA from 0.25g of soil using Qiagen DNeasy PowerSoil Kits, with a slight modification to the manufacturer's protocol (extending centrifugation time to 1.5 min after adding C3; Pulido-Chavez et al. 2022). For the four-year samples, DNA was extracted with Qiagen PowerSoil Pro kits (PowerSoil Pro, Hilden, Germany) due to the discontinuation of the Qiagen PowerSoil kits. To enhance DNA recovery, we incubated the soil sample overnight (700µl CD1 + 100µl ATL at 4°C) before proceeding with the manufacturer’s instructions. We used the primer pair 515F-806R to amplify the V4 region of the 16S rRNA gene for archaea and bacteria (Caporaso et al., 2011) and the primer pair ITS4-fun and 5.8S to amplify the ITS2 region for fungi (Taylor et al., 2016) and sequenced at the UCRInstitutee for Integrative Genome Biology (Illumina MiSeq 2 x 300 bp) using the previously published protocols (Pulido-Chavez et al,. 2022).

 We used qPCR to estimate bacterial and fungal abundance in our DNA extracts. The Eub338/Eub518 primers were used for bacteria (Fierer et al., 2005), and the FungiQuant-F/FungiQuant-R primers were used for fungi (Liu et al., 2012). The qPCR master mix was prepared with 1X of iTaq SYBR mix (Bio-Rad) and 0.4 μM final concentration of primers, with 1 μl of DNA and ultrapure water for a final volume of 10 μL.

Each reaction was run in triplicate in 384 well plates on BioRad CFX Opus 384 Real-Time PCR System starting at 94 ◦C for 5 min, followed by 40 cycles of a denaturing step at 94 ◦C for 20 s, primer annealing at 52 ◦C for bacteria or 50 ◦C for fungi at 30 s, and an extension step at 72 ◦C for 30 s. Standards were generated by cloning the 16S region of Escherichia coli with a plasmid length of 157 or the 18S region of Saccharomyces cerevisiae with a plasmid length of 353 into puc57 plasmid vectors, which were constructed by GENEWIZ, Inc. (NJ, USA) as previously established (Pulido-Chavez et al., 2022). Melt curves were generated and copy numbers were extracted using the following equation,  Gene copy number = (Starting quantity of DNA (SQ)(ng) * 6.022 * 10^{23) / (Length of target genes (bp) * 660 * 1 * 10}9) (Whelan et al. 2003), where SQ is the average of 3 technical replicates of the staring concentration (ng) value generated with CFX Maestro software and  6.022 × 10²³ is the Avogadro’s constant. We also consider the conversi^{n factor of 10}9 and the average mass of 1bp dsDN^A which is 660. 

Sequence data was processed with Qiime2 version 2020.8 (Bolyen et al., 2019). Primers were removed using cutadapt (Martin, 2011) before chimeric and low-quality regions were filtered out to produce Amplicon Sequence Variants (ASVs; DADA2 version 2020.8) (Bolyen et al., 2019). Taxonomy was assigned to each ASV using SILVA version 132 for bacteria (Yilmaz et al., 2014),  and UNITE version 8.2 for fungi (Abarenkov et al., 2020). We excluded bacterial sequences not assigned to Kingdom Bacteria and those assigned to mitochondria and chloroplasts, and fungal sequences that were not assigned to Kingdom Fungi. We assigned fungal ASVs to functional guilds using the FUNGuild database (Nguyen et al. 2016) and retained only ASVs with a “highly probable” confidence score as classified.

Statistical analyses

We used linear mixed effect models (lme4 package in R; Bates et al., 2015) with ash depth as the predictor variable to test if burn severity affected initial soil POC and MAOC concentrations 17 days after burning. The models included MAOC or POC as dependent variables and plot and soil series as random effects to account for sampling design and potential differences between soil classifications. We plotted model residuals to ensure that they were normally distributed and ran simulation-based dispersion tests to assess if data were over-dispersed (Hartig & Lohse, 2021); the linear models on untransformed variables passed both of these checks. We used the “anova” function to calculate F values, Satterthwaite’s degrees of freedom, and statistical significance (P < 0.05; “lmerTest” package in R; Kuznetsova et al. 2017).

We used a similar linear mixed effect modeling approach to assess how various dependent variables (POC, MAOC, PyOM-C, bacterial abundance, fungal abundance, bacterial richness, fungal richness, relative abundance of fungal guilds, extracellular enzyme activity, MAOC and POC δ13C values, and leaf litter decomposition rates) differed between burned and unburned plots over time after the fire. Models included fire, time, and the interaction between fire and time as fixed predictor variables. Plot and soil series were included as random effects to account for the sampling design. We also included a unique variable for each sampling location as a random effect to account for repeated measures. We plotted model residuals to ensure that they were normally distributed and ran simulation-based dispersion tests to assess if data were over-dispersed (Hartig & Lohse, 2021); potential activity of BG, NAG, and phenol extracellular enzymes was log transformed to meet assumptions of normality. The “anova” function was used to assess which predictor variables were statistically significant, as described above (Kuznetsova et al. 2017).

The effects of fire and time on bacterial and fungal community composition were assessed using permutational analysis of variance (PERMANOVA) with the “adonis” function (vegan package in R; Oksanen et al., 2013). A blocking term was included using the “strata” argument to account for repeatedly measuring the same locations. The PERMANOVA was run on a dissimilarity matrix created using the “avgdist” function in R (Oksanen et al., 2013). The matrix used the Bray-Curtis dissimilarity index and reported median distances after 100 iterations. We also used this dissimilarity matrix to conduct principal coordinates analysis (PCoA; wcmdscale function; Oksanen et al., 2013). Scores from the first and second PCoA axes were plotted to visualize community composition, and scores from the first PCoA axis were included in models predicting POC and MAOC concentration (see below).

 We used a model selection method to determine how microbial activity affected soil POC and MAOC in burned plots. Linear mixed-effect models contained potential activity of all eight extracellular enzymes as predictor variables and either POC or MAOC as response variables. We also ran separate models with fungal abundance, fungal richness, bacterial abundance, bacterial richness, fungal community composition (PCoA axis one), and bacterial community composition (PCoA axis one) as predictor variables. We did not run similar models including plant cover as predictor variables because we did not measure plant cover at the same times when soils were collected for analysis. All models included sampling plot, soil series, and sampling location as random effects to account for study design and repeated measures. To decide which predictor variables were significant, we compared all possible models after dropping each fixed term with a likelihood-ratio test (drop1 function in R) (R Core Team, 2023; Tredennick et al., 2021).