Fire alters plant microbiome assembly patterns: integrating the plant and soil microbial response to disturbance.

It is increasingly evident that the plant microbiome is a strong determinant of plant health. While the ability to manipulate the microbiome in plants and ecosystems recovering from disturbance may be useful, our understanding of the plant microbiome in regenerating plant communities is currently limited. Using 16S ribosomal RNA (rRNA) gene and internal transcribed spacer (ITS) region amplicon sequencing, we characterized the leaf, stem, fine root, rhizome, and rhizosphere microbiome of < 1-yr-old aspen saplings and the associated bulk soil after a recent high-intensity prescribed fire across a burn severity gradient. Consistent with previous studies, we found that soil microbiomes are responsive to fire. We extend these findings by showing that certain plant tissue microbiomes also change in response to fire. Differences in soil microbiome compositions could be attributed to soil chemical characteristics, but, generally, plant tissue microbiomes were not related to plant tissue elemental concentrations. Using source tracking modeling, we also show that fire influences the relative dominance of microbial inoculum and the vertical inheritance of the sapling microbiome from the parent tree. Overall, our results demonstrate how fire impacts plant microbiome assembly, diversity, and composition and highlights potential for further research towards increasing plant fitness and ecosystem recovery after fire events.

Introduction

As wildfires are currently growing in size, severity, and frequency (Westerling et al., 2006; Miller et al., 2009; Adams, 2013), understanding how vegetation regenerates post‐fire is essential in mitigating the ecological effects of fire and may have applications toward accelerating ecosystem recovery. A growing body of literature recognizes that microorganisms living inside or in close association with plant tissues, collectively known as the plant microbiome, are integral to plant health and survival (Compant et al., 2005; Santoyo et al., 2016). For instance, the plant microbiome can confer resistance to pests (Hubbard et al., 2019), provide nutrients (Moyes et al., 2016), and promote plant growth (Santoyo et al., 2016). However, an understanding of plant microbiome assembly in a post‐fire context is currently lacking. Such an understanding of the plant microbiome response to fire may be particularly important for quaking aspen (Populus tremuloides), as it is widely distributed throughout North America (Little, 1976) and aspen forests support habitat for numerous wildlife and understory plant species (Perala, 1990). Aspen is a fire‐associated species and relies largely on fire to regenerate (Shinneman et al., 2013). This regeneration occurs primarily through vegetative reproduction, specifically through root‐associated suckering, from living roots that survive even intense fires belowground (Brown & DeByle, 1987). Understanding the aspen microbiome, especially after fire, could lead to management strategies that promote the sustainability of this important plant species, as has been done in native prairie restoration (Koziol et al., 2018).

As with other Populus species (Hacquard & Schadt, 2015) the aspen microbiome and factors influencing the composition of the microbiome have been partially described. While aspen genotypes differ substantially in their ability to produce phenolic glycosides and condensed tannins (Holeski et al., 2012), which act as defense compounds against insect herbivory and fungal pathogenesis (Hwang & Lindroth, 1997; Holeski et al., 2009), the bacterial leaf endospheric community seems to be consistent across genotypes (Mason et al., 2014). However, certain bacterial strains in the aspen leaf endosphere have the ability to metabolize these compounds, possibly leading to reduced plant herbivore defense (Mason et al., 2016). In this way, leaf chemistry, the microbiome, and health of aspen are inextricably linked. Additionally, fire may alter aspen leaf chemistry, because aspen saplings growing in post‐fire environments have greater concentrations of phenolic glycosides in their leaves which may reduce herbivory (Lindroth & Clair, 2013). However, it is currently unclear whether this response interacts with the plant microbiome and if this has consequences for plant health and ecosystem recovery.

Given that fire significantly changes the soil microbiome (Xiang et al., 2014; Glassman et al., 2015; Whitman et al., 2019) and that the soil microbiome is a major source of plant endospheric microorganisms living within tissues both belowground and aboveground (Grady et al., 2019), it seems likely that plants regenerating in burned soils will have different microbiomes than those growing in unburned soils. Similarly, fire can alter numerous edaphic properties (Certini, 2005; Dove et al., 2020b), which have also been shown to impact plant microbiomes (Shakya et al., 2013). However, in the context of clonal plants, such as aspen, the plant microbiome in regenerating ramets may be somewhat resistant to the effects of fire if the source of microbial inoculum is the buried, living roots and rhizomes which are insulated from the effects of fire. It is these lateral roots that give rise to clonal saplings (or ‘suckers’) that are characteristic of aspen regeneration immediately following fires. Understanding the sources of the endospheric microorganisms in such clonal plants is fundamental to resolving mechanisms of plant microbiome assembly and the functional role of lateral roots in this process. However, the relative importance of horizontal (i.e. soil to sapling) or vertical (i.e. rhizome to sapling) transmission of microbes in determining the endospheric microbiome composition in clonal plants is currently unknown. From a biotechnology standpoint, understanding microbial assembly patterns could inform the design of successful microbial inoculation treatments, especially those intended to increase the success of plant regeneration post‐fire (Pizarro‐Tobías et al., 2015).

The current understanding of the impacts of fire on the plant microbiome is confined to the plant rhizosphere, and predominantly consists of the studies of mycorrhizal symbionts rather than the broader microbiome (Jonsson et al., 1999; Stendell et al., 1999; Glassman et al., 2015; Dove & Hart, 2017). From the mycorrhizal literature, it is becoming clear that the impact of fire on mycorrhizal colonization of plant roots is generally negative in terms of both diversity and colonization rate. For example, a recent global‐scale meta‐analysis shows that, on average, fire reduces mycorrhizal colonization by 21% (Dove & Hart, 2017), which likely negatively impacts plant survivability post‐fire. However, fire has also been shown to select for certain plant growth‐promoting bacteria in the holm oak (Quercus ilex) rhizosphere (Fernández‐González et al., 2017). Metagenomic assessment of these same samples showed an enrichment of nitrogen (N) cycling functional genes, which may further support plant growth by increasing nutrient availability (Cobo‐Díaz et al., 2015). These studies indicate that fire impacts plant rhizospheres, which may lead to changes in plant health and recovery. However, we know of no study that has investigated the response of the endospheric plant microbiome to fire and how the collection of plant microbiomes across tissue types interact with each other and the soil in a post‐fire context. Given the importance of the plant microbiome for plant health, especially in stressful conditions such as after fire, understanding the plant microbiome response to fire could constrain predictions of plant propagule success and the trajectory of plant community recovery in post‐fire environments.

To improve our understanding of the aspen microbiome responses to disturbance by fire, we sampled < 1 year old saplings or ‘suckers’ (clonal ramets) from burned and unburned areas along a burn severity gradient and sequenced microbial DNA extracted from leaves, stems, fine roots, rhizomes, rhizospheres, and bulk soils. We hypothesized that not only would the microbial community differ among habitats (e.g. leaf, stem, root, rhizome, rhizosphere, and bulk soil), but the impact of fire on the microbial community composition would vary among habitats (i.e. a habitat × fire interaction on microbial community composition). Specifically, we expected that because the plant microbiomes were not directly impacted by fire (saplings emerged in the weeks after the fire disturbance), endospheric microbiomes would be relatively less impacted by fire when compared to the rhizosphere and bulk soil microbiomes. We hypothesized that these changes in microbiome composition would be partially explained by plant and soil chemistry, which may lead to altered plant health through differences in the functional roles of these microbiomes (e.g. potential plant pathogens and mycorrhizal fungi). Using source tracking modeling (Shenhav et al., 2019), we also hypothesized that due to decreased microbial abundance in soil post‐fire (Dove et al., 2020b), the relative importance of vertical transmission (i.e. sourced from the rhizome rather than the soil) of endospheric microorganisms in fine roots, stems, and leaves would increase with burn severity. Our overall goal is to characterize the aspen microbiome and its response to fire and to elucidate how the microbiome responds to the effects of fire. As this is a novel area of research, our findings provide fundamental knowledge of plant microbiome community assembly and how this is impacted by disturbance that could be beneficial towards utilizing microbial communities to enhance plant and ecosystem regeneration.

Materials and Methods

Site description, experimental design, and sample collection

The study was conducted on the Fishlake National Forest, in central Utah, USA. Dominant tree species in our study area include trembling aspen (Populus trembuloides), subalpine fir (Abies lasiocarpa), Douglas‐fir (Pseudotsuga menziesii), and Engelmann spruce (Picea engelmannii). Soils are in the Sessions‐Faim‐Embargo‐Elwood‐Clayburn and Bickmore families and classified as Agricryolls (USDA‐NRCS, 2015). Our study area lies between 2900 and 3000 m elevation, and the mean annual temperature and precipitation for the site is 4°C and 762 mm, respectively.

As part of the Fire and Smoke Model Evaluation Experiment (FASMEE) project (Prichard et al., 2019), on 20 June 2019, Fishlake National Forest personnel conducted a high‐intensity prescribed fire that reached the forest crown (essentially a planned ‘wildfire’) across 2200 acres with the objective of studying fire and smoke behavior to improve fire‐modeling efforts, as well as the forest management objectives of reducing encroaching conifers and restoring aspen clones via increased post‐fire suckering (i.e. vertical growth from buried rhizomes), thereby improving wildlife habitat. We collected samples later that growing, season about 12 weeks after the fire, on 12–13 September 2019 when regenerating aspen were about 0.5 to 1 m in height. We sampled plants and soils from nine, 5 m radius plots in each burn severity level (unburned, moderate‐severity burn, and high‐severity burn) for a total of 27 plots (Supporting Information Fig. S1) over the three sites. These sites were relatively flat (i.e. little slope or aspect influence), on the same soil type (see earlier), and within a relatively small geographic area (< 1.5 km). At the high‐severity burn site and the unburned site, we laid out plots in two bisecting 100 m north–south and east–west transects in the shape of an ‘X’ and the sampled 5 m diameter plots spaced every 25 m across these transects. At the medium severity site, a contiguous location for fully bisecting transects could not be located due to the patchy nature of the burn in that area, so instead of an ‘X’, our transects formed a ‘T’ where only the northern end of the north–south transect intersected the middle of the east–west transect. We assessed soil burn severity visually in the field following the scale of Parsons et al. (2010; Fig. S1) and photographed representative plots. At each plot we sampled and composited plant tissues (leaves, stems, fine roots, and rhizomes) from three < 1 year old aspen saplings within a 5 m radius of the plot centroids (plant tissue types were kept separate). Rhizosphere was operationally defined and collected as the soil adhering to the fine roots at the time of collection, thus the roots were not cleaned before placing samples in sterile Whirl‐Pak (Madison, WI, USA) bags in the field. We also sampled and composited c.30 g of bulk soil near each plant (within 30 cm) to a depth of 5 cm. Plant tissues and soils were sampled into Whirl‐Pak bags, sealed, then stored and shipped on dry ice, and were kept frozen in ultra‐cold freezers (−80°C) in the laboratory until further analysis.

Sample pre‐processing and DNA extraction

Before DNA extraction, stems, fine roots, and rhizomes were washed and surface‐sterilized as previously described (see Cregger et al., 2018). Leaves were not surface‐sterilized because the freezing process in the field caused the leaves to become brittle and break, so the leaf microbiome represented a combination of leaf phyllosphere (surface) and internal leaf endosphere (tissue) microbiomes. Plant tissues were then cut into fine pieces (c. 5 mm) and extracted with the Qiagen PowerPlant Pro DNA Kit (Qiagen, Venlo, the Netherlands) following the standard protocol with the following exceptions: (1) before extraction, frozen stem, fine root, and rhizome tissues were bead‐beaten on a Retsch tissue lyser (Qiagen, Venlo, the Netherlands) for two 1‐min intervals at 30 s−1 with a sterile 3 mm steel bead, (2) 50 mg of leaf or bead‐beaten stem, fine root, or rhizome tissue samples were lysed in sterile bead tubes using the Precellys 24 tissue homogenizer (at 5500 rpm for three cycles of 30 s bead‐beat, 30 s rest; Bertin Instruments, Montigny‐le‐Bretonneux, France), and (3) 50 μl of the optional phenolic separating solution from the Qiagen PowerPlant Pro DNA Kit was added to leaf extractions during the homogenization step. Rhizosphere soil was collected as the pre‐sterilized rinsate of the fine roots which had adhering soil attached from time of collection. Rinsates were centrifuged at 10 000 , and we removed the supernatant. We then used the Qiagen PowerSoil DNA Kit (Qiagen) to extract these rhizospheres (i.e. the soil pelleted after centrifugation of the rinsate) as well as the bulk soils following the standard protocol again using a Precellys tissue homogenizer to bead‐beat extractions. We used a Zymo DNA Clean and Concentrator‐5 kit (Zymo Research Corp., Irvine, CA, USA) to purify and concentrate all extractions before polymerase chain reaction (PCR) amplification. Extractions were quantified using the Qubit dsDNA BR Assay Kit (Invitrogen, Waltham, MA, USA).

PCR amplification, sequencing, and bioinformatics

A two‐step PCR approach was used with barcode tagged‐templates and primers targeting the V4 region of the 16S ribosomal RNA (rRNA) gene for Archaea and Bacteria and the internal transcribed spacer 2 (ITS2) region for Fungi using pooled primer sets to increase coverage of archaeal, bacterial, and fungal taxa (Cregger et al., 2018 – see Supporting Information Table S1). The first step of PCR included 2.5 μM of peptide nucleotide acid (PNA) blockers targeting the plant plastid and mitochondria rRNA gene for 16S rRNA gene amplifications (GGCAAGTCTTCTTCGGA and GGCTCAACCCTGGACAG) and 2.5 μM of PNA targeting plant nuclear rRNA genes for ITS2 region (CGAGGGCACGTCTGCCTGG) were used to reduce amplification of plant material. Each reaction contained 2μl of template DNA, 0.25 μM of primer pair, 1× of KAPA HiFi HotStart ReadyMix, and molecular grade water for a total reaction volume of 25 μl. PCR amplifications were performed with the conditions 95°C for 3 min, 25 cycles (30 cycles for endosphere) of 95°C for 30 s, 78°C for 30 s, 55°C for 30 s and 72°C for 30 s and a final extension of 72°C for 5 min. The second step of PCRs were amplified following the Illumina 16S Metagenomic Sequencing Library Preparation note with the conditions 95°C for 3 min, eight cycles of 95°C for 30 s, 55°C for 30 s and 72°C for 30 s and a final extension of 72°C for 5 min.

After PCRs, all experimental units were pooled based on band intensity and purified with Agencourt AMPure XP beads (0.7 : 1 bead to DNA ratio; Beckman Coulter Inc., Pasadena, CA, USA). Paired end sequencing (2 × 251) was completed on pooled prepared libraries on an Illumina MiSeq instrument (San Diego, CA, USA) at Oak Ridge National Laboratory using V2 chemistry and included a ≥ 15% PhiX sequencing control library.

Both 16S and ITS2 datasets were denoised, joined, delineated into amplicon sequence variants (ASVs), and assigned taxonomy in the qiime2 environment (v.2019.7, Bolyen et al., 2019). Before ASV delineation using dada2 (Callahan et al., 2016), 16S reads were truncated to 200 bases (to remove low quality base calls) with the first 19 bases trimmed (to remove primers). For ITS2, reads were trimmed (including primers) using the itsxpress plugin under the default parameters (Rivers et al., 2018) with no further trimming/truncation prior to ASV delineation. We then assigned representative sequences a taxonomic classification using the Naïve Bayes classifier through the sklearn python package for 16S rRNA sequences with the SILVA database (release 132; Quast et al., 2013) and a confidence of 0.7. We assigned taxonomic classifications of the ITS2 of the ribosomal operon to representative sequences using consensus blast (percent identity: 80%; e‐value: 0.001; minimum fraction of assignments: 0.51; Camacho et al., 2009) and the UNITE reference database (v.8.0; Abarenkov et al., 2010). We removed remaining 16S reads assigned as mitochondria and chloroplasts and kept only reads assigned to Bacteria and Archaea. On average, these removed reads accounted for 67% and 10% of 16S reads in the plant tissue and soil samples, respectively. All ITS reads were assigned to the fungal kingdom. Fungal ASVs were further classified as mycorrhizal or potentially pathogenic using funguild (Nguyen et al., 2016) because the direct impact of these two guilds on plant fitness is reasonably defined. For this classification, ‘Possible’, ‘Probable’, and ‘Highly Probable’ guild classifications were used, and in instances where multiple guild classifications were given (i.e. saprotrophic/ectomycorrhizal), guilds besides mycorrhizal and pathogenic were dropped (there were no instances of ASVs having both mycorrhizal and pathogenic classifications). The funguild classifications were manually curated to remove classifications disproven in the Populus system, of which we did not find any. All ASVs assigned to the family Glomeraceae were classified as arbuscular mycorrhizal (AM) fungi.

Soil and plant chemistry

To determine relationships between microbial community composition and soil and plant chemistry, frozen sub‐samples were sent to the University of Georgia Extension Soil, Plant, and Water Laboratory for chemical analyses. Soils were dried, ground, and analyzed for total carbon (C) and N concentrations by direct combustion using the Elementar vario MAX CNS Element Analyzer (Elementar, Langenselbold, Germany). Additionally, fresh soil was measured for pH in a well‐mixed 1 : 1,w/v soil : calcium chloride (CaCl2) slurry (0.01 M) using a Fisherbrand accuTupH Rugged Double Junction pH Combination Electrode (Waltham, MA, USA). Bulk soil was additionally assessed for inorganic N concentrations by extracting 5 g of fresh soil with a 20 ml 1.0 M potassium chloride (KCl) solution. Extracts were filtered and measured for ammonium and nitrate concentrations using an Astoria continuous flow analyzer (Astoria‐Pacific, Clackamas, OR, USA) and a OI Analytical FS 3100 (College Station, TX, USA), respectively. We did not measure inorganic N on rhizosphere soils because we were limited by sample mass.

Plant C and N were assessed by direct combustion similar to soils earlier. Concentrations of boron (B), calcium (Ca), iron (Fe), potassium (K), magnesium (Mg), manganese (Mn), phosphorus (P), sulfur (S), and zinc (Zn) were determined by microwave digestion and inductively coupled plasma‐optical emission spectrometry (ICP‐OES). Briefly, dried samples were digested using 0.5 g of sample and 10 ml of nitric acid (HNO3) in a CEM Mars 6 Microwave (Matthews, NC, USA; Kingston & Walter, 1995). Digests were then measured using ICP‐OES using a Spectro Arcos FHS16 (Kleve, Germany; Creed et al., 1994).

Statistical analyses

All statistical analyses were conducted in r v.4.0.2 (R Development Core Team, 2008) with the betareg (Cribari‐Neto & Zeileis, 2010), emmeans (Lenth, 2020), hillr (Li, 2018), nlme (Pinheiro et al., 2017), phyloseq (McMurdie & Holmes, 2013), and vegan (Oksanen et al., 2013) packages. The R code used to conduct statistical analyses and generate figures can be found at https://github.com/nicholascdove/burned_aspen_microbiome.

Differences in α‐diversity were compared by means of Hill numbers (Jost, 2006) of samples rarefied to 2000 reads for 16S and 1000 reads for ITS (the average integer of reads was used after 999 rarefactions) at orders of q = 0, q = 1 and q = 2 (rarefaction curves are presented in Fig. S2). Hill numbers express the effective diversity of a sample (i.e. the number of equally abundant species that would be needed to give the same value of a diversity measure) among different metrics of q (Chao et al., 2014). Because the parameter q determines the relative weighting of rare species, multiple traditional α‐diversity indices (e.g. richness, Shannon’s, Simpson’s) can be compared in a unified framework by adjusting the q metric. For instance, at q = 0, all species are weighted equally (richness); at q = 1, species are weighted proportionally to their relative abundance (analogous to Shannon’s index); and at q = 2, rare species are down‐weighted (analogous to Simpson’s index). Differences in means of Hill numbers among habitats and burn severity were assessed by nested ANOVA with plot as random effects. Because we were primarily interested in an effect of burn severity, where we found significant interactions between habitats and burn severity, we performed individual ANOVAs and corrected the P‐values using the Benjamini and Hochberg false discovery rate adjustment (Benjamini & Hochberg, 1995). The resulting ANOVAs did not include plot as a random effect because the resulting models had one composite sample per plot (i.e. only one habitat was tested at a time). Where independent variables were significant, we assessed multiple comparisons by Tukey’s HSD (honestly significant difference) test. We used Q–Q plots and scale‐location plots to inspect normality and homoscedasticity, respectively.

Differences in the community composition of the archaeal and bacterial and fungal microbiomes among habitats and burn severity levels were assessed by permutational multivariate analysis of variance (PERMANOVA; Anderson, 2001). For the PERMANOVAs, we used Bray–Curtis dissimilarity applied to proportionally normalized data (i.e. not rarefied). Similar to our approach for α‐diversity, where we found significant interactions, we performed individual PERMANOVAs for each habitat (e.g. soil or tissue type) and corrected the P‐values as described earlier. Heterogeneity of multivariate dispersions were tested for by using the ‘betadisper’ function in vegan (Anderson et al., 2006; Oksanen et al., 2013). When broken up by habitat, the data were homoscedastic. Archaeal and bacterial and fungal community compositions were visualized using principal coordinates analysis using Bray–Curtis dissimilarity applied to proportionally normalized data.

We determined differences in plant and soil chemistry among fire severities for each habitat with ANOVA using beta‐regression to better fit the distributions of the dependent variables. Beta‐regression fits models when the data are beta‐distributed and do not pass assumptions of parametric models. Where significant, multiple comparisons among burn severities were assessed by Tukey’s HSD test. The relationship between the microbiome composition and the associated plant habitat tissue chemistry were assessed individually for each habitat using distance‐based redundancy analysis (dbRDA) with Bray–Curtis dissimilarities applied to proportionally normalized microbiome data (Legendre & Anderson, 1999). Marginal effect of each term included in the model (i.e. the significance of each term) was assessed using ANOVA.

Differences in fungal potential pathogen and mycorrhizal relative abundance and richness among habitats and burn severities were assessed by nested ANOVA with plot as a random effect following the same statistical approach as for α‐diversity. Differences in ectomycorrhizal composition among habitats and burn severities were analyzed by PERMANOVA following the same approach as for the archaeal and bacterial and fungal communities.

Finally, we used feast (Shenhav et al., 2019) to determine the relative dominance of vertical (i.e. sourced from the bulk soil and rhizosphere) versus horizontal (i.e. sourced from the rhizome) microbial transmission into the leaf, stem, and fine root microbiomes, controlling for plot. feast is a Bayesian modeling approach that estimates proportions of microbial sources in a given community by leveraging its structure and measuring the respective similarities between a sink community and potential source environments using Gibbs sampling. feast is similar, but is a faster and more accurate algorithm, than the popular microbial Source Tracker algorithm (Knights et al., 2011). Differences in the relative contributions of ‘sources’ in relation to burn severity were also assessed by individual ANOVAs for each ‘sink’ microbiome.

Results

Sequencing results

After quality and taxonomic filtering (i.e. removal of plant nuclear and organelle DNA), we sequenced 2.49 × 106 16S reads across 153 samples (nine samples were removed from statistical analyses due to low read depth primarily in endophytic samples (< 2000)), with a minimum read depth of 2221 and a maximum of 160 207. For ITS, we sequenced 8.68 × 106 reads across 161 samples (one sample was removed from statistical analyses due to low read depths (< 1000)) with a minimum read depth of 1026 and a maximum of 218 783. While differences in read depth by burn severity for individual habitats occurred (Kruskal–Wallis H test: P < 0.05, Table S2), such differences were idiosyncratic and did not explain our results, suggesting that variable read depths did not affect our findings.

α‐Diversity

Archaea and Bacteria

Both habitat and burn severity significantly affected archaeal and bacterial α‐diversity (ANOVA: P < 0.05, Table S3; Fig. 1). Across the unburned samples, habitat was a significant predictor of archaeal and bacterial α‐diversity (ANOVA: P < 0.001, Table S4), and the soil and rhizosphere were about two to six times as diverse as the plant endosphere microbiomes (see Table S5 for complete statistics).

Fig. 1

Mean (and standard error) of Hill numbers (q = 0, analogous to richness) of Archaea and Bacteria and Fungi across habitats and levels of burn severity. Letters represent significant differences among burn severity levels within a habitat (P < 0.05). Note different axis scales.

We also analyzed the differences in archaeal and bacterial α‐diversity for each habitat individually. At q = 0 (analogous to species richness), burn severity was only a significant predictor of leaf archaeal and bacterial α‐diversity (leaf: F 2,25 = 8.1, P = 0.002; all others: P > 0.05; Table S6; Fig. 1), with a 40% decrease in diversity in burned relative to unburned plots on average. At q = 1 (analogous to Shannon’s diversity) and q = 2 (analogous to inverse Simpson’s diversity), we failed to detect differences in leaf archaeal and bacterial α‐diversity among burn severity levels (P > 0.05; Fig. S3; Table S6). Hence, the decrease in archaeal and bacterial α‐diversity at q = 0 in the leaf microbiome was primarily due to a loss of rare taxa in the burned plots. Instead, at q = 1 and q = 2, burn severity was only a significant predictor of the rhizosphere and bulk soil microbiomes (Fig. S3; Table S6).

Fungi

Habitat and burn severity interacted in their effect on fungal α‐diversity (ANOVA: P < 0.10, Table S3; Fig. 1). Across the unburned samples, habitat was a significant predictor of fungal α‐diversity (ANOVA: P < 0.001, Table S4), and soil and rhizosphere microbiomes were 3–5 times more diverse than the plant endosphere microbiomes (Table S5).

Because of the interaction between habitat and burn severity, we also analyzed the differences in fungal α‐diversity for each habitat individually. At q = 0, burn severity was only a significant predictor of fungal α‐diversity in belowground habitats, with the strongest differences between the high‐severity and unburned plots (Table S6; Fig. 1). Alternatively, at q = 1 and q = 2, burn severity was only a significant predictor of fungal α‐diversity in the leaf microbiome (Fig. S3; Table S6).

Microbial community composition

Habitat explained 26.1% of the variation in the Archaea and Bacteria community composition (PERMANOVA: P < 0.001, R 2 = 0.261; Fig. S4A). Because there was an interaction between habitat and burn severity level (P < 0.001), we analyzed each habitat separately. When separated by habitat, the leaf, rhizome, rhizosphere, and bulk soil microbiomes exhibited differences among burn severity levels (PERMANOVA: rhizome – P = 0.073, others – P < 0.001; Fig. 2; Table S7). For the leaf, rhizosphere, and bulk soil, these differences were noted at the phylum and class level with greater relative abundances of Firmicutes (particularly the order Bacillales) and Gammaproteobacteria (particularly family Burkholderiaceae) in fire‐affected samples (Figs S5–S7). For the rhizome, however, taxonomic differences were noted at the genus level with greater relative abundances of Pseudomonas in burned samples and Streptomyces in the high‐severity burn samples (Fig. S8). Overall, Archaea represented relatively few reads in our dataset (0.44% of 16S reads in the final ASV table) and decreased in relative abundance in the rhizosphere and bulk soil microbiomes with increasing burn severity (Fig. S5). Among habitats, we did not detect evidence of spatial autocorrelation of the archaeal and bacterial microbiome composition at each burn site, suggesting that plots within sites were indeed independent (Fig. S9; Table S8).

Fig. 2

Principal coordinate analysis (PCoA) ordinations of Archaea and Bacteria community composition (using Bray–Curtis dissimilarity) across habitats and levels of burn severity. The percentage in parentheses represents the variance explained by each axis. The asterisk (*) next to each habitat title indicates significant differences among burn severities (PERMANOVA: P < 0.05; see Supporting Information Table S7 for full statistics). Note different axis scales.

Habitat explained 16.3% of the variation in fungal community composition (PERMANOVA: P < 0.001, R 2 = 0.163, Fig. S4B). Because there was an interaction between habitat and burn severity level (P < 0.001), we analyzed each habitat separately. Similar to patterns of the archaeal and bacterial communities, among habitats, the fungal community composition displayed differences among burn severities in only the leaf, rhizosphere, and bulk soil (PERMANOVA: P < 0.05; Fig. 3; Table S7). These differences were noted at multiple taxonomic levels. At the phylum level, increasing burn severity led to an increase in the ratio of Basidiomycota to Ascomycota reads in both the leaves and bulk soil as well as an increase in the relative abundances of Mortierellomycota and Mucoromycota in the rhizosphere (Fig. S10). At the family level, the effect of fire corresponded to an increase in the relative abundances of the plant pathogen containing families Ustilginaceae and Pleosporaceae (particularly genus Alternaria) in the leaf (Figs S11, S12). At the genus level, increased burn severity was associated with a shift in the bulk soil ectomycorrhizal community from Cortinarius‐, Inocybe‐, and Wilcoxina‐dominated to Russula‐dominated; however, these genus‐level differences in the bulk soil were only somewhat apparent in the rhizosphere (Fig. S12). Among habitats, we only found spatial autocorrelation for the fungal rhizosphere community composition in the unburned site (we did not detect spatial autocorrelation in other habitats and sites, Fig. S13; Table S8).

Fig. 3

Principal coordinate analysis (PCoA) ordinations of fungal community composition (using Bray–Curtis dissimilarity) across habitats and levels of burn severity. The percentage in parentheses represents the variance explained by each axis. The asterisk (*) next to each habitat title indicates significant differences among burn severities (PERMANOVA: P < 0.05; see Supporting Information Table S7 for full statistics). Note different axis scales.

Relationships between microbiome composition and soil and plant chemistry

Soil pH was significantly greater in the high‐severity burn treatment relative to the unburned site in both the rhizosphere and the bulk soil (Table S9). Fire did not significantly affect bulk soil C and N concentrations or these elemental concentrations in the rhizosphere (Table S9). However, inorganic N concentrations increased with increasing burn severity in the bulk soil (Table S9). In the bulk soil, differences in soil chemistry related with both the archaeal and bacterial and fungal community composition, explaining 33.5% and 27.2% of the composition, respectively (Table S10; Fig. 4a,b). For bulk soil Archaea and Bacteria, pH, total N, and ammonium cation (NH4 +) were the only significant individual terms in the dbRDA (Table S11). In contrast, for bulk soil Fungi, all terms but pH were significant in the dbRDA (Table S11). Similar to the bulk soils, in the rhizosphere, differences in soil chemistry related with archaeal and bacterial and fungal community composition, explaining 21.1% and 17.0% of the composition, respectively (statistics in Table S10; Fig. 4c,d). Soil pH, C and N were all significant individual terms for both the archaeal and bacterial (Table S11) and fungal rhizosphere communities (Table S11).

Fig. 4

Distance‐based redundancy analysis (dbRDA) ordinations of archaeal and bacterial bulk soil (a), fungal bulk soil (b), archaeal and bacterial rhizosphere (c), fungal rhizosphere (d) community composition across burn severities (using Bray–Curtis dissimilarity). Arrows represent the direction and magnitude (arrow length) of the relationship between the microbiome composition and soil chemistry. The percentage in parentheses represents the variance explained by each axis. Note: different axis scales. Key: carbon = C (%), nitrogen = N (%), ammonium = NH4 (mg kg−1), nitrate = NO3 (mg kg−1).

Among the different plant tissues, numerous elemental concentrations varied across burn severities (Table S12). However, we did not detect generalizable patterns across plant tissues. Nevertheless, similar to differences in microbiome composition, leaf tissues had the greatest differences in elemental composition with burn severity. For instance, N, P and Mn concentrations increased with burn severity while K and Ca concentrations decreased with burn severity (Table S12). Unlike the soil microbiomes, the plant tissue microbiome compositions generally did not relate to the chemistry of the plant organ (Table S10). In the singular exception, the fungal leaf microbiome, 50.2% of the variation in the leaf fungal microbiome composition was significantly explained by the elemental composition of the leaf tissue (Fig. S14). However, when assessed individually, none of the elemental concentrations were significant predictors of the leaf fungal microbiome composition (P > 0.10).

Potential pathogens and mycorrhizal fungi

Our fungal guild classifications resulted in 55 arbuscular mycrorrhizal (AM), 554 ectomycorrhizal (EM), and 243 potential pathogen ASVs (Table S13). The relative abundance of potential fungal pathogens differed among habitats and was highest in aboveground plant tissues (ANOVA: F 5,119 = 27.3, P < 0.001, Fig. S15). However, the effect of burn severity on the relative abundance of potential fungal pathogens interacted with habitat (F 10,160 = 3.6, P < 0.001), so we investigated each habitat individually. We found that burn severity only influenced the relative abundance of potential fungal pathogens in the leaf microbiome (leaf: F 10,160 = 7.7, P = 0.016; all others P < 0.05; Figs 5, S15). In the leaf microbiome, the relative abundance of potential fungal pathogens increased almost four‐fold in high‐severity plots compared to unburned plots. This increase was mainly attributed to the genera Cladosporium and Alternaria, which individually more than doubled and increased 30‐fold in relative abundance, respectively (Fig. 5). Additionally, Erysiphe was not present in the unburned leaf habitat, but made up almost 10% of the leaf microbiome under high burn severity (Fig. 5).

Fig. 5

Mean relative abundance of potential fungal pathogen reads as a proportion of all fungal reads in the leaf across levels of burn severity. Different letters represent significantly different total pathogen relative abundances among levels of burn severity (α = 0.05). Error bars represent standard error of the mean of total pathogen relative abundance. The three fungal pathogen genera with the highest relative abundances are highlighted in color.

Fire resulted in compositional changes to the EM community (PERMANOVA: P = 0.003, R 2 = 0.053, Fig. S16) without affecting the relative abundance or α‐diversity of EM reads across habitats (ANOVA: P < 0.05; Figs S17, S18). Interestingly, the effect of fire on the relative abundance of specific genera differed among habitats. For example, while fire resulted in an increase in Wilcoxina and Tuber in the fine root endosphere, these genera decreased in the bulk soil microbiomes in fire‐affect plots. In contrast, fire decreased the relative abundance of Russula in the fine roots but increased Russula abundance in the bulk soil microbiome (Fig. S16). Fire decreased relative abundance of AM reads in the bulk soil, but not in the rhizosphere or fine root tissues (Fig. S17). There were only 10 AM ASVs after rarefaction; this prevented robust diversity analyses.

Endospheric microbiome source tracking

Across the sink microbiomes (e.g. leaf, stem, and fine root) we were able to determine the source for 20%, on average, of the archaeal and bacterial and fungal communities (Fig. 6). The relative contribution of these sources was in some instances impacted by level of burn severity. The relative contribution of these sources varied with burn conditions for the archaeal and bacterial stem microbiome (F 4,47 = 2.6, P = 0.053), fungal leaf microbiome (F 4,80 = 5.6, P < 0.001), and fungal stem microbiome (F 4,80 = 5.5, P < 0.001). The contribution of rhizome‐derived Archaea and Bacteria to the stem microbiome increased (Fig. 6). Furthermore, while in unburned plots the rhizosphere was a substantial source of leaf and stem fungal taxa (~40%), the contribution of the rhizosphere to these microbiomes in burned conditions was significantly smaller (~4%, Fig. 6). The source contributions to the archaeal and bacterial leaf and fine root microbiome as well as the fungal fine root microbiome did not differ among burn severities (Fig. 6).

Fig. 6

Percentage of microbial taxa sourced from the rhizosphere, bulk soil, or rhizome across endospheric habitats and levels of burn severity for Archaea and Bacteria and Fungi. Error bars represent the standard error of the mean. Different letters represent significantly different groups (α = 0.05). Note: y‐axis varies, and percentage does not add to 100% because some microbial taxa were attributed to ‘unknown’ sources. Unknown sources represented, on average, 80% of source taxa for both Archaea and Bacteria and Fungi.

Discussion

This work demonstrates, for the first time to our knowledge, that fire impacts the broader plant microbiome, outside of the bulk soil and rhizosphere (for fire impacts on the rhizosphere see: Cobo‐Díaz et al., 2015; Fernández‐González et al., 2017). Similar to previous studies we show that the soil microbiome is responsive to fire (Xiang et al., 2014; Taş et al., 2014; Whitman et al., 2019; Yang et al., 2020); however, we build upon this by showing that certain plant tissue microbiomes also appear to change in response to fire, affecting both α‐diversity and community composition. While the plant microbiome composition in certain plant habitats was unresponsive to fire, fire explained the greatest proportion of both the archaeal and bacterial and fungal composition in the leaf phyllosphere. The impact of fire on the leaf microbiome was even greater than that on the rhizosphere or bulk soil microbiomes. This was unexpected because unlike these two soil microbiomes, the leaves, of course, emerged after the fire, so their microbiomes did not experience the direct impacts of the fire disturbance (e.g. heat‐induced mortality). These results, therefore, highlight the importance of the indirect impacts of fire on the plant microbiome as has been demonstrated in soils (Mikita‐Barbato et al., 2015; Whitman et al., 2019; Adkins et al., 2020).

Contrary to our hypothesis, the indirect impacts of fire on the plant microbiome composition were difficult to relate to specific chemical changes within the plant tissues. For example, while half of the analyzed elements differed significantly in concentrations among burn severities in the leaf, these differences only related with the fungal community composition, not the archaeal and bacterial composition. The lack of a relationship between the archaeal and bacterial composition and plant elemental concentrations was unexpected because leaf N concentration correlates positively with chlorophyll content (Niinemets, 1997). Greater chlorophyll contents and photosynthesis should provide additional sugars to support microbial growth in the leaf (van der Wal & Leveau, 2011). Indeed, in other plant species such as Diplotaxis tenuifolia, Pueraria montana, Quercus macrocarpa, Spinacia oleracea, and Tsuga spp. the leaf phyllosphere microbial community structure has been shown to be correlated with leaf N content (Jumpponen & Jones, 2010; Darlison et al., 2019; Dove et al., 2020a; Shahrtash & Brown, 2020). It is possible that post‐fire differences in the microbial community might be explained by differences in plant metabolites or site factors that were not measured, such as solar irradiation. For example, aspen produces numerous defense metabolites that may impact microbial community composition (Flores & Hubbes, 1979; Lindroth & Hwang, 1996), and production of these metabolites may be altered during high growth rate periods in exposed areas such as after fire (Donaldson et al., 2006). Specifically, phenolic glycosides in leaf tissues have been shown to increase with burn severity and increased levels of light, which may combat potential mammalian herbivory (Lindroth & Clair, 2013). The production of these defense compounds may also affect the leaf microbiome composition, as was recently shown in the rhizosphere communities of Populus trichocarpa by Veach et al. (2019). Furthermore, solar irradiation may be a better correlate of photosynthesis and sugar production than leaf N content. Future research into the plant microbiome response to fire should prioritize detailed measurements of the plant metabolome response as well as differences in site characteristics such as solar irradiation and temperature between burned and unburned areas. As the indirect effects of fire on microbial communities are generally considered to be longer lasting than the direct impacts (at least in soils, Hart et al., 2005), understanding these associations between the microbial communities and the indirect effects of fire may be useful in predicting and promoting plant sustainability post‐fire.

Unlike the plant tissue microbiomes, changes in the rhizosphere and bulk soil microbiomes with burn severity, were consistently explained by differences in soil chemistry associated with fire. Increased pH and inorganic N in our burned sites are corroborated by numerous studies with samples from vastly different ecosystems (Covington & Sackett, 1992; Neary et al., 1999; Dove et al., 2020b). Consistent increases in soil pH after fire are primarily due to the denaturation of organic acids (Certini, 2005) and the release of base cations during the incomplete combustion of organic matter (Arocena & Opio, 2003). Given that pH is consistently a strong determinant of soil microbial community composition (Fierer & Jackson, 2006; Tedersoo et al., 2014), it was not surprising that fire‐enhanced pH correlated with differences in microbial community composition among burn severities.

Nitrogen availability within a few years after fire is generally high (Covington & Sackett, 1992; DeLuca & Sala, 2006), originating from fire‐induced N mineralization (St John & Rundel, 1976), increased rates of organic matter decomposition (Kaye & Hart, 1998), and increases in N‐fixing plant abundance (Johnson et al., 2005). Increases in N availability can decrease N limitation of soil microorganisms and increase nitrification (Kurth et al., 2014; Hanan et al., 2016). Therefore, increases in inorganic N in the bulk soil should favor copiotrophic microorganisms with high N‐demand (Ramirez et al., 2012) as well as those that rely on ammonia for energy production (i.e. ammonia oxidizers). Indeed, the relative abundance of Proteobacteria, which are generally regarded as copiotrophs, increased in burned plots. However, prominent ammonia‐oxidizing bacteria (e.g. Nitrosospira and Nitrosomonas) and archaea (e.g. Nitrocosmicus and Nitrososphaera) had similar relative abundances across burn severities in our study. This could reflect the heat sensitivity of these populations (Dunn et al., 1985), which leads to longer recovery times (Yeager et al., 2005; Dove et al., 2020b). Nevertheless, our results contribute to the growing literature showing that the soil microbiome is not only sensitive to fire but also depends on the severity of the fire, in part due to changes in soil chemistry (Whitman et al., 2019).

Increased relative abundance of pathogenetic fungi, specifically Alternaria, Cladosporium, and Erysiphe, in the leaf phyllosphere with increasing burn severity could impact plant health, survivability and re‐establishment of beneficial ecosystem properties after fire. Alternaria spp. are common pathogens of aspen plants, and leaves infected with Alternaria can become discolored and eventually senesce (Dey & Debata, 2000). To combat Alternaria, aspen generate phytoalexins to suppress the germination of fungal spores and reduce infection (Flores & Hubbes, 1979). It is possible that after disturbance, increased light to the understory resulted in increased growth and photosynthetic capacity (evidenced by increased leaf N) at the expense of pathogen protection. Such trade‐offs have been shown in aspen (Donaldson et al., 2006) and are likely to impact the survivability of these aspen ramets. We were unable to confidently identify bacterial ASVs that were potentially pathogenic to aspen in our study from the existing literature. However, it is likely that such Bacteria exist in our dataset, and future pathogen screening of common Bacteria that change in response to fire will be important in understanding aspen fitness post‐fire.

Our hypothesis that fire would impact mycorrhizal associations of the aspen plants was somewhat supported by the data. Interestingly, although fire is generally considered to reduce both EM colonization and diversity (Dove & Hart, 2017), we found that the relative abundance and richness of EM reads in the root endosphere, rhizosphere, and bulk soil did not depend on burn severity. However, changes in the composition of the EM community are corroborated by numerous other post‐fire studies (Jonsson et al., 1999; Stendell et al., 1999; Glassman et al., 2015), which suggests that EM species have variable abilities to withstand the effects of fire and colonize post‐fire environments. For instance, laboratory experiments have shown that spores of Wilcoxina mikolae are relatively resistant to temperatures up to 65°C (Peay et al., 2009). Interestingly, in our study we found increased Wilcoxina relative abundances in burned fine root samples, which may be the result of heat‐stress tolerance. It is likely that trade‐offs exist between stress resistance and nutrient uptake (Grime, 1977), which could affect plant nutrient use and the degree to which plants benefit from mycorrhizal associations (i.e. along the mutualism‐parasitism continuum, Johnson et al., 1997). These trade‐offs may have manifested in increased abundance of Russula (and Basidiomycetes, in general), which typically decrease in response to fire (Pérez‐Izquierdo et al., 2021). It is possible that at our study site reduced EM competition post‐fire may have enabled proliferation of Russula, as antagonistic relationships among other EM fungi are common (Koide et al., 2005). These differences among studies highlights the nuanced response to fire of the soil microbiome among ecosystems.

It is also possible that differences in modes of microbial colonization among burn severities affected changes in the microbial community composition. Our hypothesis that vertical transmission (i.e. microbes sourced from the rhizome) would increase in dominance relative to horizontal transmission (i.e. microbes sourced from the surrounding soil) with increasing burn severity was somewhat supported by the data. While we did not measure microbial biomass, it is well‐documented that fire results in decreased microbial abundances in upper soil depths (Dooley & Treseder, 2012; Pressler et al., 2019). It is likely that this resulted in decreased inoculum for horizontal colonization. However, our source tracking modeling approach comes with several limitations. For instance, we were unable to pinpoint the source of a majority of sink taxa in any burn severity × habitat × amplicon combination. This might have occurred because sequencing depths were not sufficient to classify the source of rare taxa. Also, other sources besides the rhizome and soil exist. For example, dust‐derived microorganisms may play an important role in colonizing environments where microbial abundances are relatively low (e.g. recently burned areas; Barberán et al., 2015). Additionally, our analysis may have missed ephemeral microorganisms that were once but are no longer present in the soils post‐fire and pyro‐aeolian inoculation sources (i.e. microbes traveling on ash particles), which may play a particularly important role in post‐fire microbial colonization (Kobziar et al., 2018). Nevertheless, these source modeling results present an opportunity to develop hypotheses of how clonal plant species accumulate endospheric microorganisms. Future studies of source tracking modeling with robust spatial and temporal sampling coupled with empirical confirmation could provide invaluable information on plant microbiome assembly and on best practices for microbial inoculation.

No previous study has characterized the whole plant microbiome response to fire, and thus, our findings provide fundamental knowledge of plant microbiome community assembly and how this is impacted by disturbance. Overall, our results demonstrate a strong impact of fire on the aspen microbiome, which differs between microbial habitats in the plant and soil. It is likely that these microbiome differences contribute to the variable survival of aspen recruits post‐fire (Smith et al., 2011; Long & Mock, 2012) as differences in leaf pathogen load and EM composition are key determinants of plant health (Smith & Read, 2008; Dean et al., 2012). Furthermore, unknown interactions between plant‐associated bacteria and the plant may also influence plant survival. Future research should prioritize elucidating these microbiome effects through inoculation/sterilization glasshouse growth trials (e.g. Hewitt et al., 2016). As fires are increasing in size, severity, and frequency worldwide (Adams, 2013), our understanding of post‐fire revegetation will become increasingly important in maintaining the critical ecosystem services provided by plants.

Author contributions

NCD and CWS designed the study and collected samples. Laboratory work was conducted by NCD, DMK and AAC. NCD analyzed the data and wrote the manuscript with critical input from DMK, AAC, MAC and CWS. All authors contributed to the article and approved the submitted version.

Fig. S1 Overview of study design.

Fig. S2 Rarefaction curves across habitats and burn severity levels.

Fig. S3 Hill numbers at q = 1 and q = 2 across habitats and levels of burn severity.

Fig. S4 Principal coordinate analysis ordinations of microbial community composition across habitats and burn severity.

Fig. S5 Relative abundance of archaea and major bacterial taxonomic groups across habitats and levels of burn severity.

Fig. S6 Relative abundance of top 10 dominant bacterial orders across habitats and levels of burn severity.

Fig. S7 Relative abundance of top 10 dominant bacterial families across habitats and levels of burn severity.

Fig. S8 Relative abundance of top 10 dominant bacterial families across habitats and levels of burn severity.

Fig. S9 Spatial autocorrelation of 16S community composition dissimilarity at the site scale among different sites and habitats.

Fig. S10 Relative abundance of fungal phyla across habitats and levels of burn severity.

Fig. S11 Relative abundance of top 10 dominant fungal families across habitats and levels of burn severity.

Fig. S12 Relative abundance of top 10 dominant fungal genera across habitats and levels of burn severity.

Fig. S13 Spatial autocorrelation of ITS community composition dissimilarity at the site scale among different sites and habitats.

Fig. S14 Distance‐based redundancy analysis ordination of fungal leaf community composition across burn severities.

Fig. S15 Relative abundance of fungal pathogen reads as a proportion of all fungal reads across habitats and levels of burn severity.

Fig. S16 Relative abundance of dominant ectomycorrhizal fungal genera as a proportion of all classified ectomycorrhizal fungal reads across levels of burn severity for the fine root, rhizosphere, and bulk soil microbiomes.

Fig. S17 Relative abundance of arbuscular mycorrhizal and ectomycorrhizal fungal reads across levels of burn severity for the fine root, rhizosphere, and bulk soil microbiomes.

Fig. S18 Hill numbers of ectomycorrhizal fungi after rarefaction across habitats and levels of burn severity.

Table S1 Primer and Illumina adapter sequences for polymerase chain reaction amplification.

Table S2 Differences in average read depth among burn severities for each habitat.

Table S3 Statistics for two‐way ANOVAs of α‐diversity in response to differences in burn severity, habitat, and their interaction.

Table S4 Statistics for ANOVAs of α‐diversity in response to differences in habitat in only unburned samples.

Table S5 Multiple comparisons of differences in α‐diversity across unburned habitats.

Table S6 Tukey’s honestly significant difference test among levels of burn severity.

Table S7 Statistics for PERMANOVAs of the 16S and ITS communities among different habitats with burn severity as the independent variable.

Table S8 Spatial autocorrelation of microbial community composition dissimilarities among habitats at the site‐scale.

Table S9 Differences in soil pH, carbon (C) and nitrogen (N) across burn severities for the bulk soil and rhizosphere.

Table S10 Statistics of distance‐based redundancy analysis across amplicons and habitats.

Table S11 The F‐values for individual terms of rhizosphere and bulk soil distance‐based redundancy analyses assessed marginally.

Table S12 Plant elemental chemistry.

Click here for additional data file.

Table S13 Fungal amplicon sequence variants classified as ectomycorrhizal, arbuscular mycorrhizal, and potential pathogens.

Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.

Click here for additional data file.