DOI: 10.5061/dryad.83bk3j9pj
title: Watershed and fire severity are stronger determinants of soil chemistry and microbiomes than within-watershed woody encroachment in a tallgrass prairie system

description: 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.
description: <strong>Study
 site</strong> 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<sub>4</sub> grasses: big bluestem (<em>Andropogon
 gerardii</em>), Indian grass (<em>Sorghastrum
 nutans</em>), little bluestem (<em>Schizachyrium
 scoparium</em>) and switchgrass (<em>Panicum
 virgatum</em>). 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). <strong>Experimental
 design</strong> 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 (<em>Cornus
 drummondii</em>) and staghorn sumac (<em>Rhus
 glabra</em>) with intermittent wild plum (<em>Prunus
 americana</em>) and some other less frequent species (e.g., eastern
 red cedar – <em>Juniperus virginiana</em>; 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 <em>C. drummondii</em> 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<sup>3</sup> of
 supplemental small-diameter timber (native <em>Quercus
 </em>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
 <em>C. drummondii </em>shrub islands; ten plots assigned to
 low severity treatments within the established shrub islands; ten plots
 assigned to high severity treatments outside the shrub<em>
 </em>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. <strong>Soil temperature measurements
 before, during, and after burning</strong> 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. <strong>Soil
 sampling</strong> 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. <strong>Soil
 chemistry analyses</strong> 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<sub>4</sub><sup>+</sup>] and nitrate
 [NO<sub>3</sub><sup>–</sup>]), 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<sub>4</sub><sup>+</sup> and
 NO<sub>3</sub><sup>–</sup>) 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). <strong>DNA
 extraction, PCR amplification, and sequencing</strong>
 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<sub>2</sub>, 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<sub>2</sub>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). <strong>Sequence Data
 Processing</strong> 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 &gt;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 Naï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 (&lt;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<sub>Obs</sub>) 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.
 <strong>Statistical
 Analyses</strong> 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
 <strong>watershed</strong> and <strong>time
 </strong>(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 <strong>vegetation</strong> and
 <strong>time </strong>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 <strong>vegetation, severity,
 </strong>and <strong>time </strong>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
 &gt; 0.59; P &lt; 0.02) and heteroscedastic (Welch’s tests:
 F<sub>7,34</sub> &gt; 2.79; P &lt; 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<sup>2</sup>) 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. <strong>Research location
 description:</strong> 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 <em>Cornus
 </em>patches of between 5-10m in diameter which contained two, 1x2m
 plots. Outside shrub island / grass – refers to area
 near, but outside <em>Cornus</em> 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<sub>4</sub><sup>+</sup>] and nitrate
 [NO<sub>3</sub><sup>–</sup>]) 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.
description: 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.