Effects of short-interval reburns in the boreal forest on soil bacterial communities compared to long-interval reburns.

Increasing fire frequency in some biomes is leading to fires burning in close succession, triggering rapid vegetation change and altering soil properties. We studied the effects of short-interval (SI) reburns on soil bacterial communities of the boreal forest of northwestern Canada using paired sites (n = 44). Both sites in each pair had burned in a recent fire; one site had burned within the previous 20 years before the recent fire (SI reburn) and the other had not. Paired sites were closely matched in prefire ecosite characteristics, prefire tree species composition, and stand structure. We hypothesized that there would be a significant effect of short vs. long fire-free intervals on community composition and that richness would not be consistently different between paired sites. We found that Blastococcus sp. was consistently enriched in SI reburns, indicating its role as a strongly 'pyrophilous' bacterium. Caballeronia sordidicola was consistently depleted in SI reburns. The depletion of this endophytic diazotroph raises questions about whether this is contributing to-or merely reflects-poor conifer seedling recolonization post-fire at SI reburns. While SI reburns had no significant effect on richness, dissimilarity between short- and long-interval pairs was significantly correlated with difference in soil pH, and there were small significant changes in overall community composition.

Introduction

The boreal zone is one the world's largest biomes, spanning 1.89 billion ha across the northern hemisphere (Brandt et al. 2013). This zone consists of forests of cold-tolerant tree species, lakes, rivers, wetlands, and naturally treeless areas, such as shrublands and grasslands (Brandt 2009).

The Canadian boreal forest represents 28%—about 552 million ha—of the world's boreal zone (Brandt et al. 2013) and provides habitat for thousands of species, supplies numerous ecosystem services including timber, timber products, and water filtration, is home to 12% of Canada's population, and offers many other economic and cultural resources (Bogdanski 2008). Furthermore, this system stores about 10%–30% of the global terrestrial carbon stocks, mostly belowground in peatlands and soils (Bradshaw and Warkentin 2015, Kasischke et al. 2000), which may be threatened by changing fire regimes (Ribeiro-Kumara et al. 2020, Walker et al. 2019).

Fire is a common and widespread disturbance throughout much of the western Canadian boreal zone, where the average fire-free interval (FFI; number of years that pass between a pair of fires at the same site) has been observed to range anywhere between 30 and hundreds of years (Larsen 1997, Stocks et al. 2001). Fire is a critical event for maintaining healthy boreal ecosystems by shaping vegetation composition, soil chemical properties, and animal communities (Rowe and Scotter 1973, Whitman et al. 2018). Over the past 50 years, there has been a shift in the forest fire regime for many areas of the North American boreal forest, including lengthened burn season, increased lightning ignitions, and increased area burned (Wotton and Flannigan 1993, Kelly et al. 2013, Jain et al. 2017, Veraverbeke et al. 2017, Hanes et al. 2019). These shifting disturbance regimes can have adverse effects on ecosystems and may degrade forest resilience to fire (Johnstone et al. 2016). Forest resilience, broadly, is the ability for a forest to return to predisturbance conditions, often determined by ecological memory of past states (e.g. via seed banks; Johnstone et al. 2016) and the regeneration of plant communities (Gill et al. 2017). Under historical fire regimes, forests often have self-regulatory processes that limit fire frequency—e.g. for fire to spread, there must be adequate biomass and fuel accumulation (Peterson 2002, Héon et al. 2014). However, under drought conditions and as the forest ages, these self-regulatory processes can weaken by drying out available biomass and fuels, making accumulation of vegetation less important and allowing for decreased FFIs (Parks et al. 2018). In the relatively uncommon case when young forests (< 20 years) reburn, decreased FFIs (short-interval reburns) in boreal forests can alter vegetation composition (Whitman et al. 2019b), change aboveground plant production (Johnstone and Chapin 2006), and potentially induce forest-type conversions (i.e.Picea spp.-dominated to Populus tremuloides—or Pinus banksiana-dominated; Johnstone and Chapin 2006, Gill et al. 2017). Short-interval (SI) reburns also reduce soil organic horizon thickness (Johnstone and Chapin 2006, Hoy et al. 2016), change soil chemical properties by depleting total C and N (Pellegrini and Jackson 2020, Pellegrini et al. 2020), and can potentially decrease microbial decomposition rates (Köster et al. 2014, Pellegrini et al. 2020). Furthermore, novel wildfire regimes have been found to have long-term effects on biogeochemical soil processes, decreasing mineral soil organic carbon (SOC), soil extracellular enzyme activity, and soil microbial respiration (Dove et al. 2020), which may interact with vegetation responses to fire (Knelman et al. 2015). However, it is less clear how SI reburns within boreal forests may affect soil microbial communities.

Soil microbial communities provide numerous critical ecosystem functions, including cycling carbon and nitrogen, supporting plant growth and diversity, preventing erosion, and maintaining soil structure via biofilms and fungal hyphae (Van Der Heijden et al. 2008, Saleem et al. 2019). Fire can affect the soil microbial communities directly, via heating and oxidation of the soil environment, and indirectly, by increased exposure to climatic variation and from changes to the physicochemical environment (Hart et al. 2005). Immediately postfire, microbial biomass may decrease due to direct killing of microbes or the loss of nutrient resources (Dooley and Treseder 2012, Holden and Treseder 2013, Pressler et al. 2019). While wildfires—particularly high severity fires—might be expected to decrease microbial richness, effects of fire on richness have been inconsistent, with some studies reporting a decrease in bacterial richness with increasing fire severity (Sáenz de Miera et al. 2020), but others reporting a lack of significant changes (Pressler et al. 2019, Whitman et al. 2019a). Microbial community structures may take decades to recover to previous states, often requiring plant community to reestablish first (Dooley and Treseder 2012, Ferrenberg et al. 2013, Whitman et al. 2022). The effect of fire on microbial community composition is influenced by fire severity and fire-induced changes to vegetation, moisture, pH, and soil carbon (Whitman et al. 2019a, Hart et al. 2005, Holden and Treseder 2013, Sáenz de Miera et al. 2020, Day et al. 2019).

Microbial resistance (insensitivity to disturbance) and resilience (the rate of recovery after disturbance) may inform our understanding of ecosystem resilience following fire (Shade et al. 2012). Microbes have innate traits that differ from those of some of their larger-organism counterparts (high abundances, widespread dispersal potential, comparatively rapid growth potential, and comparatively rapid evolutionary adaptations; Shade et al. 2012). The stability of microbial communities over time can be influenced by biological attributes characteristic of individuals (e.g. dormancy or phenotypic plasticity), populations (e.g. dispersal rate or adaptability), or communities (e.g. richness, evenness, or microbial interactions; Shade et al. 2012). Microbial communities may show functional resistance by transitioning to a compositionally different community, yet remaining functionally redundant (Shade et al. 2012), in that the ecosystem process rates of interest remain unchanged (Allison and Martiny 2008). SI reburns may challenge resistance or resilience of microbial communities, by changing soil properties and vegetation community composition, both of which are factors that shape microbial communities (Chandra et al. 2016, Woolet and Whitman 2020, Van Der Heijden et al. 2008, Bardgett and van der Putten 2014).

There are few studies in any ecosystems that examine how microbial communities respond to these potentially transformative SI reburns, and most of these are limited to fungi. At the community level, Oliver et al. (2015) observed that soil fungal communities in a loblolly pine forest were significantly different from unburned controls under 1–2 years fire intervals, but not meaningfully different under 6 years fire intervals, and Egidi et al. (2016) also found that fungal community composition was influenced by prescribed fire frequency in temperate grasslands. Somewhat in contrast, Hansen et al. (2019) did not find meaningful differences in fungal abundance under different prescribed fire regimes in a pine savannah. Studies of the effects of SI reburns on soil bacterial communities remain rare, particularly in the context of natural systems such as the boreal forest. Fungi and bacteria often respond differently to wildfires—in a meta-analysis, Pressler et al. (2019) conclude bacteria are more resistant to fire than fungi, and within the domain Bacteria, the response of individual taxa to fire in the boreal forest varies widely (Whitman et al. 2019a). Because bacteria play a critical role in ecosystem functioning and structure, it is important to understand how they respond to fire regime changes.

Using a paired sample study design, Whitman et al. (2019b) showed that changes in FFIs caused changes to conifer and broadleaf recruitment, soil organic horizon depth, and herbaceous vegetation cover. Here, we set out to investigate whether there are consistent soil bacterial community responses to paired SI vs. long-interval (LI) reburns across the range of sites previously studied (Whitman et al. 2019b). Specifically, we asked (1) Do bacterial communities have a different response to SI reburns compared to normal fire intervals?, (2) Do SI reburns reduce bacterial community richness?, and (3) Which bacterial taxa respond positively to short vs. long FFIs? We hypothesized that:

There would be a significant effect of SI vs. LI reburns on bacterial community composition, where the dissimilarity of bacterial communities between paired SI and LI sites would be dependent on time since last (shared) fire (TSLF) and difference in FFI, i.e. paired sites that had a longer time to recover at the sampling date (longer TSLF) would be more similar than sites that had less time to recover, and paired sites that had greater differences in FFI would be less similar due to larger prefire differences in successional trajectory

Bacterial richness would not be consistently different between paired SI and LI sites, but would be predicted by pair-specific characteristics (TSLF, difference in FFI)

Bacteria that have higher relative abundance in SI sites than LI sites would include taxa that have previously been identified as being responsive to fires and associated with increasing burn severity in this region, such as certain taxa within the genera Blastococcus, Arthrobacter, and Massilia (Whitman et al. 2019a).

Methods

Study region description

Our study area is located in the Canadian boreal forest (boreal plain, boreal shield, and taiga shield ecozones; ESWG 1995) in the Northwest Territories and Alberta around Tu Nedhé (Great Slave Lake; Fig. 1).

Figure 1.

Location and fire history of the study area. (A) Recent (1984–2016) fire perimeters (dark gray) within the boreal forest of North America (green). (B) Fire history within the study region around Tu Nedhé (Great Slave Lake). Recent wildfires (1980–2015) are colored gray, while sampled wildfires are colored by year of occurrence. Shapes represent sampling locations. The boundary of Wood Buffalo National Park (WBNP) is outlined in gray and major roads are shown in black. (C) Detail showing an example of the sampling design of paired sites with short (triangles) and long (squares) interval sites. For example, consider the triangle (SI) and square (LI) in SW corner—both were burned in the 2004 wildfire, but only the SI site was also burned in the 1998 wildfire. Figure is modified from Whitman et al. (2019b) and is used under the Creative Commons license at http://creativecommons.org/licenses/by/4.0/.

The climate of the study area is characterized by short, warm summers, and long, cold winters. Fire frequencies in this area range from 30 years to hundreds of years between stand-replacing fires (Larsen and MacDonald 1998, Stocks et al. 2001). The landscape consists of open wetlands, forested peatlands, and forested uplands, with black spruce (Picea mariana (Mill.) Britton, Sterns and Poggenb.), white spruce (P. glauca (Moench) Voss), jack pine (Pinus banksiana Lamb.), and trembling aspen (Populus tremuloides Michx.) as the dominant tree species. The soils in this region are mostly classified as Mesisols, Gleysols, and Luvisols (Soil Classification Working Group 1998). The sample sites span a wide range of soil properties, with pH values ranging from 5.13 to 8.11 in the mineral horizons and 5.29 to 8.35 in the organic horizons, total C ranging from 0.4% (mineral horizon) to 44.6% (organic horizon), and a wide range of textures (Table 1). Detailed information on the study area and field sampling methods can be found in Whitman et al. (2019b).

Table 1.

Soil properties across drainage class, short and long FFIs, and organic and mineral soil horizons. Means reported ± standard deviation. SI = short-interval, LI = long-interval, EC = electrical conductivity, OM = organic matter, LOI = loss on ignition, and CEC = cation exchange capacity.

	Upland				Wetland
	SI organic (N = 11)	LI organic (N = 11)	SI mineral (N = 15)	LI mineral (N = 15)	SI organic (N = 7)	LI organic (N = 7)	SI mineral (N = 3)	LI mineral (N = 3)
pH	6.8 ± 1.0	6.5 ± 0.9	6.8 ± 0.8	6.2 ± 0.9	6.9 ± 0.5	7.0 ± 0.6	7.3 ± 0.5	6.9 ± 1.3
EC (mS cm–1)	0.3 ± 0.5	0.2 ± 0.2	0.3 ± 0.7	0.1 ± 0.1	0.8 ± 0.5	0.7 ± 0.4	0.2 ± 0.3	0.2 ± 0.2
% clay			9.3 ± 7.9	9.1 ± 6.0			10	18
% silt			32 ± 14	24 ± 6			26	54
% sand			59 ± 18	67 ± 9			64	28
% total N	0.45 ± 0.44	0.58 ± 0.31	0.07 ± 0.06	0.14 ± 0.14	1.69 ± 0.73	1.53 ± 0.43	0.27 ± 0.32	0.55 ± 0.69
% total S	0.11 ± 0.25	0.07 ± 0.11	0.72 ± 1.93	0.01 ± 0.03	0.39 ± 0.29	0.31 ± 0.32	0.03 ± 0.04	0.07 ± 0.10
% total C	10.5 ± 8.3	15.8 ± 8.1	1.8 ± 1.1	2.3 ± 1.6	32.4 ± 14.4	35.7 ± 9.6	4.6 ± 4.8	11.1 ± 14.5
% OM by LOI	27.4 ± 15.7	31.5 ± 15.7	5.1 ± 4.1	4.7 ± 3.6	68.6 ± 23.1	80.7 ± 8.8	4.0 ± 1.4	5.4 ± 2.0
CEC (cmol kg–1)	50.8 ± 33.8	74.5 ± 33.6	11.4 ± 5.9	12.0 ± 6.7	104.0 ± 25.2	106.1 ± 25.0	15.3 ± 9.4	28.0 ± 22.0
Ca (mg kg–1)	9480 ± 6675	9466 ± 5252	6403 ± 12 612	2051 ± 2054	20 944 ± 10 317	24 116 ± 8376	5330 ± 6660	9628 ± 13 081
K (mg kg–1)	231 ± 158	335 ± 254	72 ± 80	77 ± 55	518 ± 191	532 ± 189	139 ± 37	252 ± 184
Mg (mg kg–1)	626 ± 633	837 ± 508	262 ± 246	262 ± 224	2376 ± 1327	2262 ± 1112	344 ± 91	672 ± 653
Na (mg kg–1)	95 ± 46	91 ± 44	15 ± 8	13 ± 5	230 ± 212	232 ± 160	14 ± 15	46 ± 67

Experimental design and site assessment

In 2016, 44 sites (22 pairs) were identified and sampled; eight of these pairs were wetlands and 17 were uplands. Sites were chosen to represent the broad range of conditions characteristic of the region. Paired sites (Table S1, Supporting Information) shared the same prefire ecosite characteristics, prefire tree species compositions, stand structure, and their most recent fire (Fig. 1C), each of which burned between 1995 and 2015 (1–21 years since last fire at time of sampling). For each pair, the site with a short FFI (SI) had a FFI before the shared fire between 4 and 17 years, while the site with a long FFI (LI) had a FFI before the shared fire between 30 and 112 years, comparable to normal FFIs for this region. Fire history maps and fire-scarred trees were used to confirm dates of recent and previous fires (Whitman et al. 2019b).

Sites were selected > 100 m from roads. At each site, a 35-m transect oriented north–south was used to collect vegetation data and soil samples. Vegetation survey methods are described in detail in Whitman et al. (2019b). Briefly, along the transect, overstory structure was assessed by sampling live and dead mature trees (greater than 1.3 m in height and greater than 3 cm diameter at breast height (DBH)) using the point-centred quarter method (Cottam et al. 1953) every 5 m, recording live/dead status, tree species, and DBH to determine stem density, basal area, and species composition. For tree seedlings, we sampled stem density, species, and status (live or dead) on the east side of the transect using a 2-m wide belt transect with variable lengths depending on height classes. Seedlings, < 0.1 m in height were sampled up to 10 m along the transect, seedlings < 0.5 m were measured up to 20 m along the transect, and seedlings ≤ 1.33 m were measured along the full 35 m transect. To characterize understory vegetation, we sampled vegetation abundance of understory species and small shrubs (≤ 0.5 m) in 1 × 1 m quadrats every 5 m along the transect and estimated the % cover of exposed organic and inorganic surface substrates. At 0, 17.5, and 35 m along the transect, soil cores (5.5 cm diameter, 13.5 cm depth) were sampled by gently extruding and separating the core into organic (O) horizons (where present, up to 10 cm) and mineral (M) horizons (where present, up to 5 cm), since organic and mineral horizons generally have fundamentally different properties (Lindahl et al. 2007), occur at different depths from the surface, and would also be expected to respond to fires differently. Wetland mineral soils underlying thin organic horizons were only present within the top 13.5 cm at three sites. The three samples from the transect were pooled by horizon at each site and mixed gently by gloved hand in a bag. From these site-level samples, subsamples were collected for microbial community analysis and stored in LifeGuard Soil Preservation solution (Qiagen, Germantown, MD) in a 5 ml tube. Tubes were kept as cold as possible while in the field, then stored frozen. The remaining soil samples were air-dried and analyzed for soil properties, as described in detail in Whitman et al. (2019b; Table 1).

DNA extraction, amplification, and sequencing

DNA extractions were performed for each sample, with a blank extraction every 24 samples, using a DNEasy PowerLyzer PowerSoil DNA extraction kit (Qiagen) following manufacturer's instructions. LifeGuard Soil Preservation solution was removed from samples by thawing the sample on ice and centrifuging for 2 min at 10 000 × g. Preservation solution was gently pipetted off to remove as much as possible, and the sample was centrifuged a second time for 30 s at 10 000 × g to remove any remaining preservation solution. The sample was remixed with a spatula before weighing for extraction. Extracted DNA was amplified in triplicate PCR reactions, targeting the 16S rRNA gene v4 region (henceforth, ‘16S’). Reactions were performed in 96-well plates with each PCR mixture containing 12.5 µl Q5 Hot-Start High-Fidelity 2X Master Mix (New England Biolabs, Ipswich, MA), 1.25 µl 10 µM 515F primer (AATGATACGGCGACCACCGAGATCTACAC-barcode-TATGGTAATT GTGTGYCAGCMGCCGCGGTAA), 1.25 µl µM 806R primer (CAAGCAGAAGACGGCAT ACGAGAT-barcode-AGTCAGCCAGCCGGACTACNVGGGTWTCTAAT; Walters et al.2015, Kozich et al. 2013), 1.25 µl 20 mg ml–1 BSA, 7.75 µl nuclease-free water, and 1 µl DNA template. Positive control (bacterial isolate DNA) and negative control (nuclease-free water) reactions were included on each plate. PCR mixtures were amplified on an Eppendorf Mastercycler nexus gradient thermocycler (Eppendorf, Hamburg, Germany) at the following conditions: 98°C for 2 min, (98°C for 10 s, 58°C for 15 s, and 72°C for 10 s) × 30 cycles, 72°C for 2 min, hold at 4°C. PCR amplification success was verified via gel electrophoresis on a 1% agarose gel. The amplicon triplicates for samples and extraction blanks were pooled and normalized using a SequalPrep Normalization Plate (96) Kit (ThermoFisher Scientific, Waltham, MA), following manufacturer's instructions. Normalized samples were pooled and library cleanup was performed using a Wizard SV Gel and PCR Clean-Up System A9282 (Promega, Madison, WI). The purified library was submitted to the UW-Madison Biotechnology Center (UW-Madison, WI) for 2 × 250 paired end (PE) Illumina MiSeq sequencing.

Sequence data processing and taxonomic assignments

The University of Wisconsin-Madison Biotechnology Center performed demultiplexing on sequences. The total read count (not including blanks) was 5117 126 sequences; the minimum per sample read count was 36 594; the maximum per sample read count was 104 756; and the mean per sample read count was 61 652; sample blanks averaged 4873 reads and did not produce any visible bands in electrophoretic gels. Forward and reverse reads were imported into a Jupyter Notebook where QIIME2 (v 2019.10; Bolyen et al. 2019) and dada2 (Callahan et al. 2016) were used to filter, learn error rates, denoise, and remove chimeras, to generate operational taxonomic units (OTUs)—specifically, amplicon sequence variants. After quality control steps, a total of 3520 358 reads were retained (not including blanks). Taxonomy was assigned using the QIIME2 scikit-learn feature classifier trained on the 515f-806r region of the 99% ID OTUs (Bokulich et al. 2018) from the Silva138 database (Quast et al. 2013). OTUs that were identified as chloroplasts, mitochondria, or not Bacteria were removed (Archaea represented a mean of 0.1% of reads). All sequences are deposited in the NCBI SRA under accession number PRJNA857804.

Statistical analyses

Analyses and plotting were done with R (v.4.1.1; R Core Team 2021) in Jupyter notebooks and RStudio, using packages ‘phyloseq’ (McMurdie and Holmes 2013), ‘dplyr’ (Wickham et al. 2019), and ‘ggplot2’ (Wickham 2016). For all ANOVAs or PERMANOVAs with multiple comparisons, we used a Benjamini–Hochberg approach (Benjamini and Hochberg 1995) to control for false positives, using a false discovery rate of 0.10. Our experiment is designed to investigate trends in SI vs. LI burns across a wide range of sites characteristic of the region. Thus, we would not necessarily expect to detect trends that are unique to only a subset of the site conditions.

To test for the effect of reburn interval on bacterial community composition, we calculated Bray–Curtis dissimilarities across all sites from relative abundances, using the vegan package in R (Oksanen et al. 2018, Bray and Curtis 1957), and used permutational multivariate analysis of variance (PERMANOVA) to determine whether SI vs. LI was a significant predictor of community composition after controlling for paired sites by including site ID as a categorical variable in the model first as a predictor. We plotted the dissimilarities using NMDS and used the envfit function in the R package ‘vegan’ to map all measured site characteristics onto the ordination. To determine what factors to control for in subsequent data analyses, we used PERMANOVA to test whether drainage class and soil horizon were also significant predictors. We then calculated Bray–Curtis dissimilarities for each site-horizon pair. After controlling for soil horizon and drainage class by including these parameters in the statistical models first as predictors (except for when testing drainage class itself), we used two ANOVAs to test whether the dissimilarities for the 25 site pairs were affected by (1) TSLF or (2) difference in FFI between paired sites. Because envfit results suggested a relationship between paired community dissimilarities and pH, we also used ANOVA to test whether the dissimilarities between paired sites were significantly correlated with absolute differences in pH between sites.

Because vegetation community compositions after fire were affected by FFI in Whitman et al. (2019b), we also examined the relationship between vegetation-related parameters and bacterial community dissimilarities in LI vs. SI sites. These parameters included the effects of absolute difference in density of small understory tree stems between paired sites (including all stems, broadleaf stems, and conifer stems), difference in % understory vegetation cover, and Bray–Curtis dissimilarities of understory vegetation community.

In order to determine whether bacterial richness was affected by SI reburns, we estimated the richness of bacterial communities by using the R package ‘breakaway’ (Willis and Bunge 2015), using a weighted linear regression model (Rocchetti et al. 2011), which accounts for variance of the observations in the regression. To test for differences in richness between paired SI and LI samples in the organic and mineral horizons, we used paired t-tests. We then calculated the % difference in richness estimates between paired SI and LI sites and used an ANOVA to test for the effect of TSLF, difference in FFI, and absolute difference in total number of live stems, controlling for soil horizon and drainage class by including these parameters in the statistical models first as predictors. We also tested whether including bacterial richness improved predictions of postfire seedling recruitment by using ANOVA, with live overstory stems, site moisture, TSLF, and interval as predictor variables for the first model (as in Whitman et al. (2019b)), and live overstory stems, site moisture, TSLF, interval, and bacterial richness in the second model.

In order to determine which specific OTUs were associated with LI vs. SI sites, we estimated the differential abundances of OTUs in LI vs. SI sites using the R package ‘corncob’ (Martin et al. 2020) with the Wald testing procedure and a false discovery rate cutoff of 0.05 (chosen to be more conservative due to the high number of OTUs tested), controlling for sample pair and soil horizon by including these parameters in the statistical models first as predictors, and analyzing wetlands and uplands separately. Using the differential abundance estimates, we calculated the estimated log2-fold change in the relative abundances of the significant OTUs in the LI vs. SI sites. To identify responding OTUs that were also present in the soils of Whitman et al. (2019a), where we studied bacterial and fungal responses to fire severity in the Canadian boreal forest at sites from the same region, we used BLAST (Camacho et al. 2009). Every interval-responsive OTU from this study matched an OTU from the previous study at 99.6%–100% ID and all were classified according to their log2-fold change responses in burned vs. unburned sites from the previous study: enriched in burned sites (positive), depleted in burned sites (negative), or no significant change (neutral).

Results

Effects of SI reburns on bacterial community composition

Drainage class (PERMANOVA, P = 0.001, Fig. 2A; Table S2, Supporting Information) and soil horizon (PERMANOVA, P = 0.001, Fig. 2A; Table S2, Supporting Information) were significant drivers of bacterial community dissimilarities in the dataset. For the remaining analyses, differences between mineral and organic soil horizons and wetlands and uplands were controlled for by including these variables in regressions as detailed in the methods. SI reburns had small but significantly different bacterial community composition relative to LI reburns (PERMANOVA, P = 0.002, R2 = 0.03, Fig. 2B; Table S3, Supporting Information). Selected parameters, including % conifer seedlings, pH, and FFI were significantly associated with the second and fourth axes of the NMDS, and are plotted as arrows (envfit, P< 0.05). Dissimilarities between bacterial communities between paired SI and LI sites were significantly positively correlated with absolute differences in pH between paired sites (P< 0.001, R2adj = 0.45, Fig. 3)—i.e. sites that have bigger differences in pH have less similar bacterial communities.

Figure 2.

Nonmetric multidimensional scaling plots of Bray–Curtis dissimilarities of the bacterial community composition (k = 4, stress = 0.08) (A) First two axes. Drainage class is plotted in different panels (wetland = right, upland = left), and soil horizon is identified by shape and color (mineral (M) = gray circles, and organic (O) = orange squares). (B) Second and fourth axes. Reburn interval status is indicated by color and shape (long = green diamonds and short = orange) triangles. Lines connect site pairs. Black lines indicate vectors for selected significant (P < 0.05) site characteristics—pH, % conifer seedlings, and FFI.

Figure 3.

Relationship between Bray–Curtis dissimilarity between microbial communities for paired SI and LI sites and absolute difference in pH between paired LI and SI sites. Gray circles represent mineral horizons and orange squares represent organic horizons. Black line represents line of best fit (y = -1.95+3.89x; R2 = 0.45; P < 0.001).

Neither of our factors of interest had significant effects on bacterial community dissimilarities between paired SI and LI sites. There was not a significant effect of TSLF on dissimilarities between bacterial communities from SI and LI sites (P = 0.19; Figure S1 and Table S4, Supporting Information), and bacterial community dissimilarities between paired SI and LI sites were not clearly affected by differences in FFI (P = 0.12; Figure S2 and Table S4, Supporting Information)

Paired sites where the LI site had more total understory stems and a higher proportion of conifer stems than the SI site had more dissimilar bacterial communities (Fig. 4; Table S4, Supporting Information). Dissimilarities between bacterial communities in paired SI and LI sites were not affected by difference in total understory vegetation cover (P = 0.91; Figure S3 and Table S4, Supporting Information) or understory vegetation community dissimilarities (R2 = -0.01, P = 0.96; Figure S4 and Table S4, Supporting Information).

Figure 4.

Bray–Curtis dissimilarities of bacterial community composition between paired SI and LI reburn sites vs. the log difference (log(LI—SI) of all understory stems per hectare, understory conifer stems per hectare, and understory broadleaf stems per hectare faceted by drainage class. Soil horizon is identified by shape (mineral (M) = circle and organic (O) = square). Note log scale on x-axis. Values at zero had no stems recorded at either site at the time of sampling, and therefore, had no differences in stem count.

SI and LI differences in richness

Bacterial richness was greater in SI sites for some pairs (up to 93% greater richness in SI sites) and lower in SI sites for others (up to 53% lower richness in SI sites), but was not consistently significantly different (Fig. 5, paired t-test, P = 0.81, P = 0.39, P = 0.66, and P = 0.96, for upland mineral soil, upland organic soil, wetland mineral soil, and wetland organic soil, respectively).

Figure 5.

Boxplots of estimated bacterial community richness from paired (joined by gray lines) short and long FFI samples in the (A) upland organic soil, (B) upland mineral soil, (C) wetland organic soil, and (D) wetland mineral soil horizons (paired t-tests, P = 0.81, 0.39. 0.66, and 0.96, respectively). Each point represents a sample.

None of our factors of interest had significant effects on bacterial community richness between paired SI and LI sites. Percent change in estimated richness between paired SI and LI samples was not associated with differences in FFI between paired sites (P = 0.08; Table S5, Supporting Information). Differences in estimated richness between paired SI and LI samples were also not associated with TSLF (P = 0.80; Table S5, Supporting Information), or difference in number of stems (P = 0.11; Table S5, Supporting Information). Bacterial richness did not improve predictions of seedling recruitment, as estimated by live understory stem density (P = 0.32; Table S6, Supporting Information). After identifying differentially abundant bacteria (results below), we also wanted to know whether including differentially abundant OTUs of interest would improve predictions of seedling recruitment. None of our top responding OTUs (identified as Blastococcus, Rhizobiaceae, and Caballeronia; discussed in the following section) significantly improved predictions of the recruitment of conifer and broadleaf seedlings in an ANOVA, after accounting for live overstory stems, site moisture, TSLF, interval, and OTU abundance as predictor variables.

Differentially abundant bacteria between SI and LI burns

We identified OTUs that were differentially abundant in LI vs. SI sites in uplands (Fig. 6A; Figure S6, Supporting Information) and wetlands (Fig. 6B; Figure S5, Supporting Information). Of the most differentially abundant OTUs, one OTU from the genus Blastococcus had 6.6 and 4.7 times greater relative abundance in SI sites than LI sites, in both uplands and wetlands, respectively, and was abundant across samples (uplands: 1.5% ± 1.6%; wetlands: 0.7% ± 1.2%). Of the OTUs that were most differentially abundant in LI sites over SI sites, for uplands, we identified an OTU that had 3.9 times higher relative abundance in LI sites and was also relatively abundant overall (0.75% ± 0.51%) from the family Burkholderiaceae. For wetlands, we identified an OTU that was 5.0 times more abundant in LI sites and was relatively abundant overall (0.30% ± 0.26%) from the family Rhizobiaceae.

Figure 6.

Log2-fold change in relative abundance in long (LI) vs. short (SI) FFIs for (A) upland and (B) wetland sites, after controlling for Site ID and soil horizon by including these parameters in the statistical models first as predictors. Each point represents a single OTU. Color and shape indicate the response the OTU had in Whitman et al. (2019a; green upward triangle = positive response to fire, yellow circle = neutral response to fire, and pink downward triangle = negative response to fire). The x-axis label indicates the finest-scale taxonomy available. Size of points is scaled by the average relative abundance of OTUs in short and LI sites. Only significantly differentially abundant OTUs (pFDR< 0.05) and responses greater than |1| are plotted. Solid line indicates no difference in relative abundance between long vs. SI sites, therefore, points above the line indicate OTUs that were significantly more abundant in LI upland or wetland sites, and points below the line indicate OTUs that were significantly more abundant in SI upland or wetland sites.

Discussion

Short- and LI communities have distinct responding bacteria

While the fire ecology of plants, and even many fungi, is well-characterized and has long been studied (Cooper 1961, Seaver 1909), the bacterial equivalents of fire responder species like P. banksiana or the fungus Pyronema are only just being established. By identifying specific bacteria associated with SI or LI sites, we sought to both expand our understanding of fundamental questions about fire ecology for bacteria, as well as to offer potential hypotheses about the effects that changes to soil bacterial communities might have on ecosystem functioning. For example, taxa that are more abundant in SI sites may be able to thrive under the conditions characteristic of frequently burned areas. Thus, we might predict these microbes would exhibit traits that allow them to consume pyrogenic organic matter more effectively, thrive in areas that are more drought-prone, or be common in areas with less vegetation cover and tree regeneration. Taxa that are more abundant in LI sites could be general fire responders that are characteristic of normal fire regimes; these taxa may thrive in postfire environments, just not under the environmental conditions resulting from SI reburns, which are significantly different from more typical postfire environments (Whitman et al. 2019b, Donato et al. 2016).

Even though our study sites spanned a wide range of postfire conditions, we were able to identify several taxa that were consistently differentially abundant in SI and LI sites. Because of this experimental design, we believe that these observed responses may be robust within this region, rather than site-specific anomalies. Our most prominent responder was a Blastococcus OTU, which was the most abundant responding OTU in both wetlands and in uplands (Fig. 6). It was over four times more abundant in SI soils than in LI soils, and consistently represented a large proportion of the total community across sites. This was the same OTU (100% identical over the sequenced region) as an OTU in Whitman et al. (2019a) that was enriched in burned sites and was increasingly abundant in sites with increasing fire severity. In holm-oak forests of Spain, members of Blastococcus also increased in abundance after wildfires in rhizosphere soils vs. neighboring unburned areas (Fernández-González et al. 2017; Cobo-Diaz et al. 2015). The sequence for this taxon was also 100% identical to a recently classified Blastococcus deserti sp. isolated from a desert sample (Yang et al. 2019). While clearly a fundamentally different ecosystem, in their study, Yang et. al. (2019) noted that this isolate was able to survive at temperatures up to 50°C and utilize D-salacin as a C source, which is a compound found in the bark of Populus and Salix species (Palo 1984), when other related Blastococcus strains in their study could not. Both of these traits offer mechanisms that could help this taxon thrive in the postfire environment: its potential ability to survive high temperatures during fires, and also metabolize compounds characteristic of the organic horizon of an aspen-dominated forest—a common successional species after SI reburns in this region (Whitman et al. 2019b, Johnstone and Chapin 2006). Together, these studies highlight traits that allow us to propose this Blastococcus OTU as a classical ‘pyrophilous’ bacterium, which thrives in burned soils and even becomes increasingly abundant with more frequent and more severe burns.

For wetland sites, the most prominent LI associated taxon was classified as belonging to the family Rhizobiaceae. This OTU was detected in the Whitman et al. (2019a) study, but did not show a negative or positive response to fire in that dataset. Members of family Rhizobiaceae are often associated with nitrogen fixation (Spaink et al. 1998), which could be an important source of N in postfire ecosystems (Turner et al. 2019), particularly after severe fires when N-fixing plants increase in abundance (Smithwick et al. 2005). Some members of the Rhizobiaceae family also perform other steps of the N cycle, including denitrification (Rich et al. 2003). While these roles can be seen to be relevant postfire, it is not immediately clear why these taxa were more enriched in the LI wetland sites. One possibility is that if the SI wetland sites were drier, then the LI sites may have had conditions more conducive to denitrification (low oxygen due to high moisture, and high organic matter; Martínez-Espinosa et al. 2021) that selected for these taxa. Similarly, N fixation could also be more likely to be supported under these conditions: it is a highly energy intensive process requiring both lots of C for energy and at least locally anoxic conditions to protect the oxygen-sensitive nitrogenase enzyme (Smercina et al. 2019). However, even though this OTU from the family Rhizobiaceae is 100% identical over the sequenced region to taxa that are known to be able to perform these functions, there may still be important functional differences. We should consider these speculations on their role in postfire N cycling at LI sites as hypotheses for future testing.

For upland sites, the most prominent LI associated taxon, identified as being from the Burkholderiaceae family, also matched an OTU in Whitman et al. (2019a) that was found to be more abundant in burned sites, to the point of not being detected at any unburned sites. One possible explanation for this association could be that members of the genus Burkholderia, from this same family, have been found to be more abundant in lower pH soils (Stopnisek et al. 2014), and we observed here that SI reburns were associated with higher pH (Fig. 3). Supporting this observation, the corresponding OTU from Whitman et al. (2019a) clustered in a module associated with fire-responsive taxa that were more abundant at lower pH sites in a co-occurrence network.

This LI-associated organism may also play an important role in postfire plant establishment. While the finest-scale matching taxonomy in the SILVA database placed this LI-associated OTU within the Burkholderiaceae family, using NCBI BLAST, we found it is a 100% ID match for a Caballeronia sordidicola that was isolated from a lichen in Svalbard Archipelago (AF512827.1; Kim et al. 2017) and for a C. sordidicola spruce tree endophyte from a sub-boreal forest in British Columbia, Canada (MG561776.1; Puri et al. 2018). Furthermore, the authors of this second study identified that this isolate has the capacity to fix nitrogen (Puri et al. 2018), readily colonized pine and spruce seedlings, increasing their biomass production 4–7-fold (Puri et al. 2020a) and provided more than 50% of both lodgepole pine (Pinus contorta var. latifolia) and white spruce seedlings’ N requirements (Puri et al. 2020b). This is strongly consistent with the observation that conifer seedlings were more abundant at LI sites. While our study design does not allow us to separate the cause and effect between the increased presence of putative endosymbiont diazotroph C. sordidicola and more abundant conifer seedlings at LI sites, it clearly raises pressing questions. For example, linking our observations, we might ask whether increases in pH due to SI reburns shift the soil environment away from conditions that are optimal for C. sordidicola, exacerbating other factors that contribute to poor conifer seedling establishment, such as depleted seed banks (Johnstone et al. 2016) or suboptimal seedbed conditions associated with SI reburns (Whitman et al. 2019b).

Shorter FFIs alter bacterial communities, reflecting changes in vegetation and soil properties

In addition to the clear differences in bacterial community composition between different drainage classes and soil horizons (Fig. 2A), SI reburns also shifted bacterial community composition significantly (Fig. 2B). Although this effect was statistically significant, it explained only a small amount of the variation in bacterial community composition. The bacterial community changes between SI and LI were associated with decreased conifer dominance at SI sites (Figs 2B and 3). This echoes the types of changes that have been seen in the aboveground communities in this region (Whitman et al. 2019b), indicating that soil bacterial community composition will also be affected by predicted changes in this region—specifically, decreasing FFIs—although the relatively small effect size may suggest that soil microbes will be less dramatically affected than plants, likely requiring further investigation. The extent to which changes in the bacterial communities are being driven by changes in the vegetation community, and vice versa, is not possible to determine with this dataset and experimental design. As discussed in more detail below, we did not find a strong association between aboveground vegetation and soil bacterial communities. This lack of clear coupling may suggest that the primary drivers of changes to plant and bacterial communities with shorter FFIs differ.

Larger changes in pH were also associated with larger shifts in community composition from LI to SI sites (Figs 2B and 3). For cases where pH increased in SI burns, this could be explained if SI reburns were either higher severity, increasing total combustion and ash production, which would result in an increase in pH, particularly for combustion at higher temperatures (Bodí et al. 2014), or just simply from the compound effects of the two most recent fires, maintaining higher ash levels between the two. For cases where pH decreased in SI burns, this could potentially also be explained if the SI burns were higher severity. First, if more of the O horizon was burned away in the SI fire than the LI fire, decreased surface roughness could result in increased erosion and loss of the ash layer, decreasing the pH effect from the most recent burn at the SI sites. The pairs that had lower pH at the SI site than the LI site tended to be those with higher pH at LI sites to begin with (somewhat tautologically). Second, if carbonates were contributing to the higher pH at both sites before the burn, and these carbonates were more volatilized during high severity burns at the SI sites, this could potentially also explain the decrease in pH at SI sites. The association with soil pH changes is notable, since soil pH is regularly found to be a powerful predictor of bacterial community composition in regional datasets (Rousk et al. 2010, Bahram et al. 2018). Given this observation, it would be interesting to determine whether pH shifts are the primary driver of shifts in bacterial community composition with shorter interval reburns. Since this was not our primary question for this study, future studies could be designed to test this hypothesis.

We found that the effects of SI reburns on paired site dissimilarities were not clearly moderated by either of our predicted factors—TSLF or difference in FFI. We had predicted that increasing TSLF would decrease Bray–Curtis dissimilarities between paired sites, as communities recover and converge on similar states post fire. If this had been the case, it might have suggested that bacterial communities are resilient to the effects of SI reburns. However, our data do not indicate that increasing TSLF allows communities to converge meaningfully, over the range of TSLF studied (1–21 years). This is consistent with other studies of postfire recovery in microbial communities, which have indicated that it can take decades to over a century for soil microbial communities to return to their preburn states (Dooley and Treseder 2012, Ferrenberg et al. 2013, Dove and Hart 2017). Thus, it may not be surprising that the effects of SI reburns did not seem to decrease community dissimilarity over the range of TSLF in this study. The lack of significant correlations with TSLF and bacterial community dissimilarities does not necessarily indicate that ecological memory is not relevant for soil microbes (e.g. in terms of material legacies in the form of surviving bacteria). Rather, it could simply suggest that recovery is so slow that the range of TSLF (between 1 and 21 years) we studied was not sufficient to detect meaningful differences in the role of ecological memory in LI vs. SI sites. Indeed, while we did not detect a significant effect, we might speculate that such an effect could still emerge over longer timescales. LI and SI pair dissimilarities at longer TSLF (12–21 years) ranged from 0.46 to 0.78 (more similar), while dissimilarities at shorter TSLF (1–2 years) ranged from 0.58 to 0.91 (less similar). While the ranges overlap, they do suggest a trajectory toward more similar communities after longer times since last fire. Future studies might be designed to target an even broader distribution of TSLF across the SI–LI pairs to directly test this question, or to trace these same sites over time. However, such a sample set may be difficult to collect, as it becomes increasingly difficult to conclusively identify FFI as we move beyond the range of modern satellite data records.

Since SI vs. LI was a weak but significant predictor of bacterial community composition (Fig. 2B), we expected that communities might become increasingly dissimilar as differences in FFI increased. As differences in FFI increase, the LI site will have had increasingly longer to recover from the previous fire, so factors that could affect microbial community composition, such as soil properties or vegetation composition, may also become increasingly different. This would increase differences in the states of the two sites’ ‘ecological memories’ at the time of shared fire and extend microbial community dissimilarities between the paired sites. However, larger differences in FFI were not associated with more dissimilar communities. This suggests that the differences in FFI do not scale consistently with differences in the direct or indirect effects of fire on microbial community composition.

Failing to detect a significant correlation with community dissimilarity between paired sites and differences in FFI does not necessarily indicate that the concept of ‘ecological memory’ is not relevant for microbial communities, but it may indicate that the timescales over which it is relevant are different. For example, if dispersal of bacteria from unburned areas (whether adjacent sites or from deeper soil horizons) occurs more readily than for plants (and possibly fungi; Gill et al. 2022), that could effectively rapidly ‘reseed’ the ecological memory at both SI and LI sites.

Our finding that FFI was not correlated with community dissimilarity between paired sites may also reflect the range of possible outcomes from a SI fire. On the one hand, high-severity SI fires can occur when a recent burn has left high fuel loads, resulting in significant loss of biomass and SOM. At the other extreme, there can also be low-severity SI fires where frequent recent fires lead to surface fires with less tree mortality and SOM combustion. This range of possible outcomes from SI fires could obscure relationships with FFI.

Bacterial richness is not consistently affected by shorter FFIs

We expected that bacterial richness would not be clearly affected by SI reburns, largely because soil bacterial communities are extremely rich to begin with and, despite advances in sequencing technologies, methods for accurate richness estimation in high-throughput amplicon sequencing datasets remain limited (Willis 2019). This expectation was supported, with richness differences ranging from 93% higher (SI richer than LI) to 53% lower (LI richer than SI). While more than half of LI sites were richer than their SI pairs (Fig. 5) this difference was not significant and the % difference in richness ranged widely. Although fires often reach temperatures at the soil surface that would be expected to kill all bacteria (Pingree and Kobziar 2019), heating attenuates rapidly with depth in the soil, meaning that lethal temperatures in the surface of the soil may be accompanied by almost no changes just a few centimeters below the surface. Thus, it may not be surprising that the effects of fire on bacterial richness have been inconsistent across studies (Sáenz de Miera et al. 2020, Pressler et al. 2019, Whitman et al. 2019a). Our findings indicate that, although SI vs. LI have different bacterial community compositions (Figs 2 and 3), this difference is likely driven more by changes in the overall community structure or by a few community members, rather than significant differences in richness. While richness has often been linked to ecosystem multifunctionality (Delgado-Baquerizo et al. 2017) or resistance and resilience (Shade et al. 2012), these relationships are not always straightforward in soil communities, which generally have high functional redundancy (Nunan et al. 2017). This functional redundancy is perhaps also reflected in our finding that including bacterial richness in our models did not improve predictions of understory stem regeneration. Thus, even for the paired sites where bacterial richness was lower in SI than LI, we should not necessarily predict functional limitations.

Future directions

There are numerous future directions for this work. First, directly investigating some of the hypotheses generated here is of interest, including the extent to which bacterial community shifts due to shorter FFIs are driven by changes in soil pH. Research is needed to elucidate the relative influence of postfire microbial communities in determining vegetation community development and facilitating or impeding tree recruitment and, conversely, further research is needed to identify the influence of postfire vegetation on soil microbial community development and the strength of these interacting effects.

While identifying specific fire-responsive taxa is a critical first step in developing a fire ecology framework for bacteria, the next steps will be to investigate the functional traits of these responding taxa in a mechanistic way, in order to further understand their individual ecology and what effects changes in their populations might have on ecosystem processes. For example, the abundant and strongly SI-associated Blastococcus OTU is a general fire-responder in this region and in other parts of the world. Future experiments could determine whether this is due to an ability to survive higher temperatures and an ability to metabolize compounds characteristic of the organic horizon of an aspen-dominated forest. A second example is the depletion of diazotrophic conifer endophyte C. sordidicola in SI reburn sites. This observation raises questions about whether the depletion of C. sordidicola is contributing to—or merely parallels—poor conifer seedling recolonization postfire at SI reburns. Together, such approaches will help us develop a more robust framework of fire ecology for microbes.

Data availability

Code used for all analyses can be found at https://github.com/WhitmanLab/WoodBuffalo2016. All sequences are deposited in the NCBI SRA under accession number PRJNA857804.

Click here for additional data file.