Fire can impact terrestrial ecosystems by changing abiotic and biotic conditions. Short fire intervals maintain grasslands and communities adapted to frequent, low-severity fires. Shrub encroachment that follows longer fire intervals accumulates fuel and can increase fire severity. This patchily distributed biomass creates mosaics of burn severities in the landscape—pyrodiversity. Afforded by a scheduled burn of a watershed protected from fires for 27 years, we investigated effects of woody encroachment and burn severity on soil chemistry and soil-inhabiting bacteria and fungi. We compared soils before and after fire within the fire-protected, shrub-encroached watershed and soils in an adjacent, annually burned and non-encroached watershed. Organic matter and nutrients accumulated in the fire-protected watershed but responded less to woody encroachment within the encroached watershed. Bioavailable nitrogen and phosphorus and fungal and bacterial communities responded to high-severity burn regardless of encroachment. Low-severity fire effects on soil nutrients differed, increased bacterial but decreased fungal diversity and effects of woody encroachment within the encroached watershed were minimal. High-severity burns in the fire-protected watershed led to a novel soil system state distinct from non-encroached and encroached soil systems. We conclude that severe fires may open grassland restoration opportunities to manipulate soil chemistry and microbial communities in shrub-encroached habitats.

Study site

The study was conducted at the Konza Prairie Biological Station (KPBS, 39˚05’ N, 96˚35’ W), which hosts a Long–Term Ecological Research (LTER) site representative of native tallgrass prairie in the Flint Hills of KS, USA. The site spans 3,487 ha and remains undisturbed by row-crop agriculture. The vascular flora at KPBS (Towne 2002) is dominated by native C₄ grasses: big bluestem (Andropogon gerardii), Indian grass (Sorghastrum nutans), little bluestem (Schizachyrium scoparium) and switchgrass (Panicum virgatum). The Flint Hills and KPBS are characterized by shallow soils overlaying chert-bearing limestones and shales (Ransom et al. 1998). Our study site is classified as typical chernozem according to the Food and Agriculture Organizations (FAO) soils classification used by the United Nations. Topographic relief divides the landscape into upland plateaus with shallow soils, slopes with outcrops of limestone, and lowlands with deeper alluvial and colluvial soils. January mean temperature is –3˚C (range –9 to 3˚C), and the July mean is 27˚C (20-33˚C). Annual precipitation averages 835mm, 75% of which falls in the growing season between April and October. Historic fire intervals in the region ranged from three to four years (Knapp et al. 1998). Current conservation and rangeland management regimes use more frequent burning to suppress woody encroachment that threatens the tallgrass prairie ecosystems (Kollmorgen and Simonett 1965; Ratajczak et al. 2014a).

Experimental design

More than 20% of the 24ha fire-protected watershed was covered by woody shrubs in the early 2000’s (Briggs et al. 2005), but woody cover has since increased and exceeded 50% at the time of burning (Ratajczak et al. 2014a). The invading woody shrubs are mainly roughleaf dogwood (Cornus drummondii) and staghorn sumac (Rhus glabra) with intermittent wild plum (Prunus americana) and some other less frequent species (e.g., eastern red cedar – Juniperus virginiana; Briggs et al. 2005). While invasion of woody species leads to greater root biomass and increased heterogeneity of soil nutrients (Pärtel et al. 2007), we did not expect shifts in ectomycorrhae composition because these shrubs form arbuscular mycorrhizae (Mino 2016) like graminoids.

Within the fire-protected watershed, we selected a total of ten established C. drummondii shrub islands (ranging from ~5m to 10m in diameter) and set up two plots (1m x 2m) within each island: one assigned to low severity fire treatment (treatment 1 in Fig. S1) and the other to high severity treatment (treatment 2). We also established similar plots outside of established shrub islands that were subjected to low (treatment 3) and high severity fire treatment (treatment 4). The high severity fire treatments were created by amending plots with 0.6m³ of supplemental small-diameter timber (native Quercus spp.), whereas low severity treatments received no additional fuel. To compare the effects of the woody encroachment to a non-encroached watershed, we selected another, adjacent 35.9ha watershed that had been annually burned each spring for the past three decades. In this non-encroached watershed, we randomly selected ten plots for comparison between the two watersheds with distinct fire histories (treatment 5). This experimental design resulted in a total of fifty experimental units: ten plots assigned to high severity treatments within the established C. drummondii shrub islands; ten plots assigned to low severity treatments within the established shrub islands; ten plots assigned to high severity treatments outside the shrub islands; ten plots assigned to low severity treatments outside the shrub islands; and ten plots in the annually burned watershed assigned to low severity treatments (Fig. S1). If not ignited by the broadcast burn, the additional fuel in the high severity treatments was directly ignited. This was usually necessary within the shrub islands because the fire rarely carried through them.

Soil temperature measurements before, during, and after burning

We randomly selected one plot from each of the five treatments to measure heat penetration during fire treatments. At each of the selected plots, we dug a 20cm deep pit to install soil thermocouples. A type K thermocouple (stainless steel sheathed 24 gauge; Omega Engineering, Inc., Stamford, Connecticut, USA) was inserted at 2-, 5-, 10-, and 20cm depths into the wall of the pit in mid-March 2017 in anticipation of the burn treatment as weather conditions would permit. Each thermocouple was attached by a PVC-insulated, type K extension wire to Omega OMPL-TC data loggers, which were placed in waterproof cases containing desiccant packs and buried in playground sand approximately 2m outside the plot edge. Soil was replaced into the pit after thermocouple installation. Data Loggers were set to record every 5min starting on March 15, 2017. Both the fire-protected watershed and the annually burned watershed were broadcast burned on April 15, 2017 as part of scheduled KPBS watershed management. Data Loggers were collected on April 20, 2017.

Soil sampling

We sampled all fifty plots four days before and 28 days after the fire. Our previous fire-response experiments indicated that responses to fire are rapid and fungal communities turn over within a few weeks following a fire event (Reazin et al. 2016). After removing litter or residual coals when present, we sampled two 15cm deep mineral soil cores from each plot with a 6.35 cm diameter slide hammer impact soil corer with 15cm plastic liners (AMS Inc., American Falls, Idaho). Samples were placed in plastic bags and transported on ice to the laboratory where they were manually homogenized and subsampled for soil chemistry and nucleic acid analyses. After homogenization, two subsamples were transferred into 50 ml Falcon tubes (Corning Inc., Corning, New York) and stored at -20°C until further processing.

Soil chemistry analyses

When adequate soil homogenate was available, a 50g subsample was analyzed for soil chemistry at the Kansas State University soil testing laboratory (www.agronomy.k-state.edu/services/soiltesting/). One of the two frozen subsamples was thawed and immediately oven dried overnight at 60°C and ground to pass through a 2mm sieve. This subsample was further divided for analyses of pH, soil organic matter (SOM), total carbon (TotC), total nitrogen (TotN), total phosphorus (TotP), readily available inorganic nitrogen (ammonium [NH₄⁺] and nitrate [NO₃^–]), and plant available phosphorus (Mehlich P). A 10g subsample was used to measure soil pH directly in a 1:1 soil slurry in deionized water. One gram of dry soil was used to estimate organic matter content through loss on ignition as described in Combs and Nathan (1998). Total C and N were measured using a LECO TruSpec CN combustion analyzer (LECO, St. Joseph, MI) on a weight percent basis from a 0.45g subsample of prepared soil. Total phosphorus was analyzed colorimetrically using a modified Kjeldahl digestion and a flow analyzer from a 1g subsample of prepared soil. Inorganic nitrogen (NH₄⁺ and NO₃^–) was colorimetrically estimated from a 2g subsample extracted with 1M KCl and cadmium reduction for nitrate (Gelderman and Beegle 1998) and run in separate channels in a flow analyzer to measure the ions simultaneously. Plant available phosphorus was estimated from a 2g subsample using Mehlich 3 soil test extractant (Mehlich 1984) as described in Frank et al. (1998).

DNA extraction, PCR amplification, and sequencing

Total genomic DNA was extracted from ~0.5g soil subsamples using E.Z.N.A Soil DNA Kit (Omega Bio-Tek, Norcross, Georgia) following the manufacturer’s instructions and stored at -20ºC until PCR amplification. The DNA was quantitated using an ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, Delaware) and standardized to 2ng/µl for PCR amplification. For fungi, we chose the Internal Transcribed Spacer 2 (ITS2) region of the ribosomal RNA gene (Schoch et al. 2012) for our analyses.  We used the fITS7 (Ihrmark et al. 2012) and ITS4 (White et al. 1990) primers with unique 12bp barcodes in each 5’-end in 50µl PCR reactions. The volumes and final concentrations of reagents were as follows: 10µl forward and reverse primer (1µM), 10µL template DNA (2ng/µl), 5µl dNTP (200µM), and 0.25µL (1/2 unit) Phusion Green Hot Start II DNA polymerase (ThermoScientific, Pittsburg USA), 10µL of Phusion 5X HF Buffer with 7.5mM MgCl₂, and 14.5µL molecular grade water. The PCR began with an initial denaturing step for 30s (98°C) and was followed by 35 cycles with 10s of denaturing (98°C); 30s of annealing (54°C); 1min of extension (72°C); and, concluding with a 10min final extension (72°C). For bacteria, we targeted the highly variable V4 region of 16S ribosomal RNA gene with forward primer 515f and reverse primer 806r (Caporaso et al. 2012) appended with 12bp barcodes in each 5’end in 50µl reactions. The reaction volumes and conditions were identical to those used for fungi except for the 55°C annealing temperature and 30 PCR cycles used to generate bacterial amplicons.

For fungi and bacteria, amplification of PCR contaminants was determined by a negative PCR control in which templates were replaced with ddH₂O. Each sample was PCR-amplified in triplicate and 30µL of each amplicon was combined into one per experimental unit. The pooled 90µl amplicons were purified using the Mag-Bind RxnPure Plus Magnetic Bead Clean-up solution (Omega Bio-Tek, Norcross, Georgia) following a modified manufacturer protocol with a 1:1 ratio of PCR product to magnetic bead solution and two rinse steps with 80% ethanol. Following cleanup, a total of 200ng of amplified DNA per experimental unit was pooled. Because the negative controls yielded little measurable DNA, the entire elution volume from the cleanup (40µl) was included. Illumina adapters and indices were added using four PCR cycles, KAPA Hyper Prep Kit (Roche, Pleasenton, California), and 0.5µg starting DNA. The library was sequenced (2 x 300 cycles) using the Illumina MiSeq Personal Sequencing System at the Integrated Genomics Facility (Kansas State University, Manhattan KS USA).

Sequence Data Processing

Sequence data were processed using the mothur pipeline (v. 1.38.0; Schloss et al. 2009) following the MiSeq standard operating protocol to generate OTUs (Kozich et al. 2013). Paired fastq files were contiged and sequences with >1bp difference with the primers, without an exact match to the sample-specific identifiers, were omitted. Bacterial and fungal sequences were analyzed independently. Bacterial sequences were aligned against a mothur-formatted 16S Silva Alignment (v. 132; www.arb-silva.de), whereas fungal sequences were truncated to the length equal to the shortest high-quality read (236bp excluding primers and sample-specific identifiers). Bacterial and fungal data were then pseudo single-linkage clustered (99%; Huse et al. 2010) to control for platform generated errors and screened for potential chimeras (UCHIME; Edgar et al. 2011), with putative chimeras culled. The sequence data were assigned to taxon affinities using mothur-embedded Naïve Bayesian Classifier (Wang et al. 2007). Bacterial sequences were screened against the Ribosomal Database Project’s 16S reference training set (v.9) and sequences without a reliable assignment to Domain Bacteria (<100% bootstrap support) or assigned to Eukarya, Archaea, mitochondria, or plastids were removed. Fungal data were assigned to taxa using the UNITE taxonomy reference (Kõljalg et al. 2013) and sequences with no match in the UNITE reference or assigned to Protista and Plantae were removed. A sequence distance matrix was generated for aligned bacterial sequences. Fungal sequences were clustered using vsearch (Rognes et al. 2016). Both bacterial and fungal data were assigned to Operational Taxonomic Units (OTUs) at 97% similarity or greater. Rare OTUs represented by fewer than ten sequences in the entire dataset were removed as potential artifacts (Brown et al. 2015; Oliver et al. 2015).

We estimated Good’s coverage (i.e., complement of the ratio between local singleton OTUs and the total sequence count) for each experimental unit to evaluate the representativeness of our sampling. To estimate richness and diversity, we iteratively (100 iterations) calculated observed (S_Obs) richness, diversity (i.e., Shannon-Weiner diversity [H’]) and evenness (i.e., Shannon’s equitability) with the data subsampled to 1,000 sequences per experimental unit, as recommended by Gihring et al. (2012) to avoid biased comparisons of diversity and richness estimators in samples with unequal sequence yields. Although five times more bacterial (5,000 sequences per sample) than fungal (1,000 sequences per sample) sequence data were subsampled, coverage was lower for bacteria (0.85 ± 0.03; mean ± std. dev.) than fungi (0.95 ± 0.02). Fungal coverage estimates were generally high indicating reasonable sampling of soil communities regardless of shallow subsampling.

Statistical Analyses

Unless otherwise noted, ANOVA was performed for each research question. To test the effects of fire history across watersheds (Q1), we compared only non-encroached plots assigned to low-severity fire treatments before and after the fire in two watersheds with distinct fire histories – one burned annually, another protected from fire for more than 25 years (treatments 3 and 5 in Fig. S1). These analyses included main effects watershed and time (i.e., before and after fire) and their interaction. To test the effect of woody encroachment (Q2), we compared low severity fire treatment plots within and outside shrub islands in the fire protected watershed before and after the fire (treatments 1 and 3). These analyses included main effects vegetation and time and their interaction. To test the effect of fire severity and its dependency on the vegetation context (Q3), we compared plots in a fire-protected watershed representing two vegetation types (shrub-encroached vs. grass-dominated) and two fire (low vs. high severity) treatment combinations (treatments 1, 2, 3, and 4). These analyses included main effects vegetation, severity, and time as well as their two- and three-way interactions.

The soil chemistry data, as well as both bacterial and fungal richness and diversity data were non-normal (Shapiro-Wilk Goodness of Fit tests: W > 0.59; P < 0.02) and heteroscedastic (Welch’s tests: F_7,34 > 2.79; P < 0.04). As such, these data were ln-transformed, except for the proportional community data and percent data (OM, total N, total C) that were arcsine-square root transformed.

To visualize and infer compositional differences within fungal and bacterial communities, we calculated pairwise Bray-Curtis distances and visualized these data with Principal Coordinates Analysis (PCoA) using R package ‘vegan’ (Oksanen et al. 2019). The optimal number of dimensions (k) was selected based on stabilizing stress less than 0.20 using 1000 runs with empirical data and a random seed starting value. Community data were compared using a nonparametric permutational analog of traditional analysis of variance (Oksanen et al. 2019) using function ‘adonis’, the ‘vegan’ implementation of PERMANOVA (Anderson 2001). To address questions about community convergence or divergence we used the function ‘betadisper’. Similarities between samples were calculated for PERMANOVA using the function ‘betadiver()’ method “beta.z”. Function ‘envfit’ was used to explore goodness of fit (r²) correlations of soil chemistry variables with PCoA axes created for pairwise distances among soil microbiome samples (Table S1). To identify OTUs associated with observed community differences, if any, we performed indicator species analysis with package ‘indicspecies’ (Cáceres and Legendre 2009) using function ‘multipatt’ to test associations between taxa and treatments for each of our three research questions. Results and significance (alpha=0.05) for the one hundred most abundant bacterial and fungal OTUs are presented (Table S2-4). To visualize trends in taxonomic level indicators of change, heatmaps were created using package ‘ampvis2’ (Andersen et al. 2018). Heatmaps (Figs. S3-5) included the ten most abundant bacterial and fungal OTUs by percent read abundance within samples corresponding to each research question.

Research location description:

The study was conducted at the Konza Prairie Biological Station (KPBS, 39˚05’ N, 96˚35’ W), which hosts an annually burned (SpB) and a fire-protected (20B) watershed.

Variable description: (each variable name and description, units, missing value code)

• Non-encroached watershed – Watershed SpB
• Encroached watershed – Watershed 20B
• Within shrub island / shrub – refers to area within Cornus patches of between 5-10m in diameter which contained two, 1x2m plots.
• Outside shrub island / grass – refers to area near, but outside Cornus patches which contained two, 1x2m plots.
• High severity fire / high burn – The high severity fire treatments were created by amending plots with 0.6m3 of supplemental small-diameter timber (native Quercus spp.).
• Low severity fire / low burn – low severity treatments received no additional fuel.
• pH of soil
• soil organic matter (SOM) % arcsine-square root transformed.
• total carbon (TotC) % arcsine-square root transformed.
• total nitrogen (TotN) % arcsine-square root transformed.
• total phosphorus (TotP) ppm ln transformed.
• readily available inorganic nitrogen (ammonium [NH₄⁺] and nitrate [NO₃^–]) ppm ln transformed.
• plant available phosphorus (Mehlich P) ppm ln transformed.
• observed (SObs) richness – we calculated iteratively (100 iterations)
• diversity (i.e., Shannon-Weiner diversity [H’]) – we calculated iteratively (100 iterations)
• evenness (i.e., Shannon’s equitability) – we calculated iteratively (100 iterations)
• Operational Taxonmic Unit (OTU) – Sequence data were processed using the mothur pipeline (v. 1.38.0; Schloss et al. 2009) following the MiSeq standard operating protocol to generate OTUs (Kozich et al. 2013). OTU counts were made into proportional data based on total OTUs within each sample and then arcsine-square root transformed.

Datasets - one is for bacteria (KonzaFire_Bacteria_1_23_20) and one for fungi (KonzaFire1_10_20). One .R file with annotations, #s, and notes on data filtering, etc. The other datasets include partitioned data that is amenable to corresponding R code and functions that require manipulation of format/organization of data (descriptions of use is described in annotations in R code).

Three samples in the fungal dataset (T032tag6030.KNZFIRE; T000tag6050.KNZFIRE; T032tag6050.KNZFIRE) were mismanaged and contaminated during analysis and are therefore missing rows in the fungal dataset.