Response of fungal communities to fire in a subtropical peatland

Wildfire, an increasing disturbance in peatlands, could dramatically change carbon stocks and reshape plant/microbial communities, with long-lasting effects on peatland functions. Soil fungi are important in controlling the belowground carbon and nutrient cycling in peatlands; however, the impact of altered fire regimes on these fungi is still unclear. We assessed fungal abundance, composition, and diversity across four soil depths (0–5 cm, 6–10 cm, 11–15 cm, 16–20 cm) under low-severity and high-severity fire in a subtropical peatland in the southeastern USA. Low-severity fire significantly increased fungal Shannon diversity and saprotrophic fungi in the 0–5 cm soil layer immediately after fire and then retracted within 2 years. This pattern was not observed below 5 cm soils. The dominant fungal class − Archaeorhizomycetes declined initially and then returned to pre-low-severity fires levels at 0–5 cm depths. Time since low-severity fire was a primary driver of fungal composition in the 0–10 cm soil depth, while spatial distance among sites affected the deeper soils (11–20 cm). The fungal Shannon diversity failed to recover in the unburned state even 30 years after high-severity fire, especially in 6–20 cm soil layers. Stratification patterns of the fungal community were diminished by high-severity fire. Soil properties (either phenolics or carbon) were the primary drivers in shaping fungal community reassembly after high-severity fire across all soil depths. Collectively, the fungal communities seem to be highly resilient to low-severity fire, but not to high-severity fire in the shrub-dominated coastal peatlands.

important natural disturbance in peatlands (Pellegrini et al. 2018;Whitman et al. 2019), and its regimes have been shifting in the face of rising temperatures and increased drought events resulting from climate change and human activities (Ali et al. 2012;Langner and Siegert 2009;Turetsky et al. 2015). For example, some low-latitude fire-resistant peatlands became fire-prone over the last decade due to high-frequency low-severity fires and occasional high-severity fires (Cole et al. 2019;Turetsky et al. 2015). Fire severity and frequency changes may have cascading consequences for long-term carbon dynamics in peatlands Pellegrini et al. 2018;Turetsky et al. 2015). Therefore, it is critical to understand the consequence of fire regimes in terms of ecosystem structure, function and dynamics (Pressler et al. 2019).
Microbes play a central role in soil carbon turnover Winsborough and Basiliko 2010), and their responses to fire may substantially influence peatland recovery after a burn (Allison and Martiny 2008). Fungi, as the predominant peat decomposers in the unsaturated upper peat, principally moderate decomposition and nutrient cycling in peatlands (Phillips et al. 2013;Wang et al. 2021). Most studies have assessed the effects of highintensity/high-severity fires on the fungal community, which cause large declines in fungal biomass and result in fungal community shifts (Dooley and Treseder 2012;Holden et al. 2013;Semenova-Nelsen et al. 2019). Although devastating peat-consuming highseverity wildfires have been increasing due to current climate change, most fire events in peatlands are low-severity (Boby et al. 2010;Turetsky et al. 2015), which actually benefit carbon accretion in some peatlands (Flanagan et al. 2020). Highseverity wildfires tend to have greater impacts and longer lasting effects than low-severity fires due to fire-induced tree mortality and initiation of smoldering fires that cause deep soil combustion. Previous studies showed that high-severity wildfire had negative effects on soil fungal diversity (Day et al. 2019;Hernandez-Rodriguez et al. 2013;Holden et al. 2013;Martin-Pinto et al. 2006). While some studies observed neutral effects of increasing fire severity on fungal community in boreal wetlands . Xiang et al. (2015) even found both low-and high-intensity fire had similar negative effects on fungal alpha diversity.
Fires also modify soil fungal communities that could take years to return to its pre-fire state, thus possibly impacting their ecological functions in peatlands, such as carbon sequestration (Clemmensen et al. 2015), ectomycorrhizal colonization of plants (Treseder et al. 2004) and peat decomposition (Holden et al. 2013). An understanding of the effect of time-since-fire on fungal community is therefore a foundation for unearthing the long-term effects of fire on soil carbon in peatlands. Severity of fire and fungal community legacies could impact fungal community fire return interval. For example, low-severity prescribed fires had a minimal or no effect on microbial community (Kranz and Whitman 2019;Oliver et al. 2015), while fungal communities took at least 24 years to return to pre-fire levels after high-severity fires (Holden et al. 2013). Resistant fungi might exhibit a neutral response to a light perturbation and maintain their biomass, abundance, and composition, while resilient communities may change a little and then rapidly return to their original structure (Allison and Martiny 2008) after low-severity fire. Different fungal species or functional groups also have various recovery timespans, ectomycorrhizal and saprobic fungi differing substantially in the time required to recover to pre-fire levels (Holden et al. 2016;Treseder et al. 2004;Yang et al. 2020). In this context, understanding the influence of variations in fire severity on fungal community are important for predicting microbial feedbacks to carbon accumulation in firefrequented ecosystems.
Heat penetration and organic matter loss varied with depth during fire, where the top horizon was more affected than lower horizons (Certini 2005). Some studies observed a depth-dependent response pattern of fungal communities to wildfire (i.e., highseverity fire), a stronger response in the top horizon than in lower horizons, in boreal and cold-temperate forests (Holden et al. 2013;Yang et al. 2020). While Kranz and Whitman (2019) reported that low-severity fire did not impact bacterial community across all horizons, but they did not assess the fungal community. The degree of soil heating and organic matter combustion also depended on the fire severity (Kranz and Whitman 2019), which would cause a divergent response of fungal communities in different soil depths as heat is attenuated to a greater or lesser degree. However, these relationships are largely unknown. A recent study did show that the response pattern of soil carbon accumulation to fire depended on depth and fire-intensity (Sawyer et al. 2018).
Due to increasing shrub and tree encroachment in boreal peatlands , understanding the response of fungi to changes in fire regimes in subtropical shrub peatlands has significant global implications. However, studies of fungal community as a function of time since fire in subtropical shrub peatlands are relatively rare. Here, we study the response of fungi to low-severity fire and high-severity fire in shrub-dominated peatland in Pocosin Lakes National Wildlife Refuge (PLNWR), North Carolina, USA ( Figure S1).
In this study we address three questions: 1) how does fire severity alter soil fungal community composition, diversity, and abundance? 2) does the recovery of the fungal community vary with fire severity? and 3) does depth-dependent fungal community patterns diverge in response to variations of fire severity? We hypothesized that fungi are highly resilient to lowseverity fire as fungi possess heat-resistant structures (Greene et al. 2010), or can survive in spore banks (Glassman et al. 2016). While fungal community may lose their resilience under high-severity fire due to high temperature and destruction of the organic matter (Barcenas-Moreno and Baath 2009). We further postulated that low-severity fires may have stronger effects on surface soil than subsoil fungal communities, while high-severity fire may reset fungal communities and diminish the fungi stratified pattern in peatlands. It has been found that low-severity fire only burn the top few centimeters of the soil surface and samples that are taken down to greater depths may dilute the signal of the response of soil microbial communities to fire (Neary et al. 1999). High-severity fire can consume entire organic horizons leaving large ash deposits (Hammill and Bradstock 2006). Nutrients derived from ash that percolate into the soil (Sawyer et al. 2018) are likely to homogenize the accessible substrate for the fungal community, which may, in turn favor a similar fungal community down soil depths.

Sampling methods
In this study, space-for-time substitution (i.e., time gradients/a fire chronosequence) was used to assess the impact of fire on fungal community succession. It is a common way to study fire impact on microbial succession (Pellegrini et al. 2018;Bowd et al. 2019;Whitman et al. 2019;Yang et al. 2020). Indeed, in many situations, it may be the only way of gaining a long-term perspective on ecological processes (Cutler et al. 2017). Although this approach has disadvantages, it allows investigation of the response of real intact ecosystems to environmental factors and disturbances with high accuracy as "time-for-time" predictions (Blois et al. 2013;Bragazza et al. 2015).

Study sites and prescribed burn
The study site is in an evergreen shrub bog community at the Pocosin Lakes National Wildlife Refuge,North Carolina,. Dominant plant of study area included Ilex glabra, Pinus taeda, Woodwardia virginica and Vaccinium formosum. The mean annual temperature is 16.8 °C, and mean annual rainfall is 1,230 mm. The region was re-classified as subtropical by Trewartha in 1968 from Koppens original 1900 classification to adjust for both the changing temperature in the last century and the thresholds separating wet and dry climates (Trewartha and Horn 1980). A prescribed fire took place on Pungo North and Pungo South (2000 × 800 m) in March 2015. There was no heavy precipitation or extreme drought before the prescribed fire day. Climatic conditions were slightly wetter than normal before the prescribed fire day, as would be expected to allow fire to propagate yet remain under control. Long-term average rainfall during the months of January through March was 29.49 cm in nearby Plymouth, NC (PLYMOUTH 5 E, NC US). The total for those same months in 2015 was 31.04 cm. The U.S. Drought Monitoring Map (https:// droug htmon itor. unl. edu/ DmData/ TimeS eries. aspx) showed normal conditions in the days before the prescribed fires. To characterize the temperature and duration of prescribed fire in peatlands, we monitored soil temperatures at the surface of litter (litter charred but not the peat soil beneath) during a prescribed fire. The prescribed fire displayed a rapid consumption of aboveground biomass by fire, and soil surface temperatures reached a peak temperature of 450 °C within 3 min and then dropping below 60 °C within ten minutes after the fire front passed 1 3 (Flanagan et al. 2020), therefore, we classified this prescribed fire event as low-severity fire. Soils were collected before fire samples (i.e. undisturbed sample, UD), 15 days post-low-severity-fire (15d) and 71 days post-low-severity-fire (71d) from Pungo North and Pungo South as low-intensity fire samples. One adjacent 2 years post-low-severity-fire (2y) plot was also collected as low-severity fire samples. The 2y site had experienced a similar fire pattern as the current prescribed fire with a brief in ten minutes light surface fire (USFWS, personal observations). Based on fire records, two other sites with more than 1-m of peat burned 30 years ago and 5 years ago were selected for high-severity fire sampling, respectively. We collected samples from those two sites (5y and 30y) as high-severity fire samples ( Figure S1). The distance between those five sites ranged from 0.25 km to 3.3 km. To remove spatial effects, all of study sites are in near proximity and have similar pocosin topographic conditions and shrub/scrub plant communities. They experience same climatic condition and peat development rates thus removing climate variations among selected sites.
Based on field measurements, 5 cm was the maximum depth of heating recorded in our prescribed fire site (Flanagan et al. 2020). Given that the impacts of low-severity fire on soil rarely penetrate deeper than a few centimeters of soil layers documented by our measurement and other studies (Kennard and Gholz 2001;Sawyer et al. 2018;Flanagan et al. 2020), only the top 20 cm of peat were used for studying the heatinduced changes in fungal community. We took triplicate soil cores at each site (with a distance > 10 m from each other), and each soil core was sliced to four sub-samples (0-5 cm, 5-10 cm, 10-15 cm, 15-20 cm). We recorded shrub diversity in each plot. A total of 108 soil samples (4 depths × 3 cores × 9 site-fire time = 108) were collected. Soils were transported to the laboratory on ice. Half of the samples were frozen at − 80 °C for DNA isolation; the rest were stored at 4 °C for chemical analysis. Climate and some basic information of each site (based on different fire history) was summarized in Table S1.

Soil chemistry analysis
Dissolved organic carbon (DOC) and total soluble phenolics were measured by adding 40 ml deionized water to 2 g fresh peat samples, which were shaken on shaker for 10 min. The supernatant was filtered through a 0.45-μm membrane filter. DOC was determined by a total C analyzer (Shimadzu 5000A, Kyoto, Japan) with a method detection limit (MDL) of 0.16 mg/L. Following the Folin-Ciocalteu procedure, soluble phenolics were measured. Inorganic nitrogens (NH 4 -N and NO 3 -N + NO 2 -N, 2 M KCl extraction) were determined colorimetrically on a flow-injection analyzer (Lachat QuikChem 8000, Milwaukee, WI, Wisconsin, USA). Total carbon and nitrogen in soil were analyzed by a combustion CN soil analyzer equipped with a TCD detector (Thermo-Quest Flash EA1112, Milan, Italy). The water content was determined by weighing of peat before and after drying to constant weight at 80 °C.

DNA Extraction and Sequencing
Genomic DNA was extracted from 0.25 g (fresh weight) of each homogenized soil sample using the PowerSoil DNA isolation kit (Mo Bio Laboratories, Carlsbad, CA, USA), according to the manufacturer's instructions.
A two-step PCR modified from Lundberg et al. (2013) was used for DNA amplification. The fungusspecific primers ITS1F and ITS4 (White et al. 1990) were used. In the first step, primer constructs were used, which included a frameshift section (six per gene region), a linker section recognized in the second step, and gene specific sequences. The reaction condition was 95 °C for 10 min, 30 cycles at 95 °C for 1 min, 50 °C for 1 min, 72 °C for 1 min, and 72 °C for 10 min. The second step ligated MID tags, incorporating a linker region recognizing that used in the first step and Illumina adaptor sequences. The reaction condition was 95 °C for 10 min, 5 cycles at 95 °C for 1 min, 63 °C for 1 min, 72 °C for 1 min 20 s, and 72 °C for 10 min. A PCR negative control was performed with sterile H 2 O as template. The amplicon concentration of each sample was determined after purification using Qubit® 2.0 Fluorometer (Invitrogen, Grand Island, NY, USA). Samples were pooled at equimolar concentrations, purified using AMPure Bead cleanup, and submitted to the core facility at Duke University (Durham, NC, USA) for sequencing using Illumina MiSeq (Illumina, San Diego, CA, USA). The sequence data of 30y sites and beforefire sites were retrieved from another of our papers (accession number SRP122579, Wang et al. 2021).

3
Sequence data obtained from this study has been deposited in the NCBI Sequence Read Archive (SRA) with the accession number SRP176379.

Quantitative PCR
Quantitative PCR (qPCR) was performed to estimate the abundances fungi in soil samples. Fungi were targeted with ITS1f (5'-TCC GTA GGT GAA CCT GCG G-3') and 5.8 s (5'-CGC TGC GTT CTT CAT CG-3') primer pair (Fierer et al. 2005). Each 25-µl reaction mixture contained 12.5 µl of SYBR™ Green PCR Master Mix (ThermoFisher Scientific, MA, USA), 9.5 µl of PCR water (Mo Bio Laboratories), 1 μl of each of the primers (10 μM), and 1 μl of isolated DNA as the template. The conditions were 15 min at 95 °C, followed by 40 cycles of 95 °C for 1 min, 53 °C for 30 s, and 72 °C for 1 min. Calibration standards and negative controls were included on each plate in triplicate. A plasmid standard containing ITS region was generated using DNA extracted from pure culture strain. Genomic DNA was sourced from a Aspergillus sp. isolate. The amplification of ITS region was cloned using the TOPO TA cloning kit (Invitrogen). Plasmids were isolated using the Qiaprep Plasmid Miniprep kit (QIAGEN, Valencia, CA, USA) with DNA concentrations determined by Qubit® 3.0 Fluorometer (Invitrogen, Grand Island, NY, USA). Calibration standard curves were generated using triplicate tenfold dilutions of plasmid DNA. We used at least three nonzero standard concentrations per assay. Target copy numbers for each reaction were calculated from the standard curves, assuming that the average molecular mass of a double-stranded DNA molecule is 660 g mol −1 .
Bioinformatics ITS sequences, including from this study and retrieved from SRP122579, were quality filtered and processed using the standard QIIME pipeline ). All sequences were demultiplexed into samples based on unique identification barcodes using the split_libraries_fastq.py script with a quality score less than 20. As ITS1f and ITS4 are too long to pairend, the sequences from both primers were processed separately. Sequences were sorted by gene region and primers were removed using cutadapt, and then reads were trimmed to 220 bp using -fastq_filter with a q10 minimum value as a quality filter. Sequences were dereplicated and singletons were removed. The ITS1 and ITS2 region from the remaining sequences was identified and extracted using the fungal ITSx software (Bengtsson-Palme et al. 2013). The extracted sequences were binned into operational taxonomic units (OTUs) at 97% match using USEARCH (Edgar 2010). OTUs were classified taxonomically using a QIIME-based wrapper of BLAST against the UNITE database (Abarenkov et al. 2010) as a reference, and a BLAST e-value < 1e was a criterion for matching.
Taxonomic-based alpha diversity was calculated as Chao1 index (richness) and Shannon's diversity index (H'). All singleton OTUs were excluded from the dataset. The results from both primers reads were similar ( Fig. 1 and Figure S2), thus libraries from ITS1f reads were used for further analysis of fungal communities. In order to analyze fungal diversity and composition at the same sequencing depth, samples were rarefied to the minimum number (680) of sequences among samples for the further analysis.
To be more accuracy, we manually checked the guilds results based on published databases (Põlme et al. 2020;Smith and Read 2008;Sun et al. 2016;Tedersoo et al. 2014;Thormann 2006;Thormann and Rice 2007) and relevant published literature. All OTUs that not match taxa in the database were classified as "unassigned".

Statistical analyses
The study was divided into two datasets that were analyzed separately. Dataset1 included 84 samples was treated as low-severity fire samples. Dataset2 included 48 samples was used to investigate the response of fungal community to high-severity fire. For each dataset, the data analysis was performed separately by soil depth.
Linear mixed models with time since fire as a continuous variable were used to determine the general relationship effects of time since fire on soil properties, fungal copy numbers, fungal diversity and relative abundance of fungal groups with site as a random 1 3 factor. Time since fire and other variables were log or square root transformed when necessary. Models were fit with the package lme4 (Bates et al. 2015). Firstly, corrected Akaike information criterion (AICc) was used to identify the best mixed-effects model from linear, quadratic and cubic polynomial models. We then calculated P, marginal R 2 (R m 2 ) and conditional R 2 (R c 2 ) for each best model by the function 'r.squaredGLMM' in MUMIN (Barton 2009). Marginal R 2 (R m 2 ) represents the variance explained by fixed effects, whereas conditional R 2 (R c 2 ) represents the variance explained by both fixed and random effects. Tukey's HSD test with time since fire as categorical variable was used to compare time effects for each time of post-fire.
Principal Coordinates Analysis (PCoA) was performed in vegan to illustrate the community shifts of fungi along fire stage and soil depths. PER-MANOVA with ADONIS function, ANOSIM (analysis of similarity) and multiple-response permutation procedure (MRPP) with 999 permutations were used to evaluate the significant difference in the community structure between fire stages and peat depths. All of these analyses were performed in R with the vegan package.
To quantify the drivers behind fungal community shifts following fire and identify the effects of time since fire on fungal composition when accounting for the spatially nested design and soil heterogeneity, generalized dissimilarity modeling (GDM) (Ferrier et al. 2007) was performed using the gdm package. Due to combustion of the complete plant aboveground biomass in our study site (Flanagan et al. 2020), shrub diversity of 15d and 71d were not surveyed. Therefore, we did not include shrub diversity for low-severity fire for GDM analysis.
To determine which fungal OTUs differed between fire stages, we used the function 'manyglm' in mvabund (Wang et al. 2012). The species that were not observed in at least 60% of the samples were removed before comparison. Sequence counts were Log10 + 1 transformed, which improved the distribution of the data towards normality. We modelled sequence counts on a negative binomial distribution in the GLMs. P-values were obtained using 999 bootstraps of residuals and adjusted for multiple comparisons using a standard step-down resampling procedure.

Soil physicochemical properties
Regardless of fire severity, soil total carbon (TC) and ammonia concentration (NH 4 -N) were not affected by time since fire in any horizon, as analyzed by linear mixed-effect models (Table S2). Phenolics in 0-5 and 6-10 cm layers exhibited a quadratic relation with time lapses since low-severity fire, namely, a decrease at 15d (log e (days since fire) = 2.7) and an increase from 71d to 2y (i.e. log e (days since fire) = 4.2 and 6.6). The pH value also had a significant non-linear relationship with time since low-severity fire across all soil layers with lowest value at 71d ( Figure S3). The response of pH to low-severity fire lagged behind that of phenolics concentration in top two layers. After high-severity fire, soil properties below 5 cm soil layers had similar fluctuation pattern as time progressed for both phenolics and pH. Both measures declined through 30y at 6-10 cm layer, but increased at 30y at the two layers below 10 cm depth (Table S2 and Figure S4).

Fungal copy number
Fungal ITS copies in all soil samples amounted to 10 8 -10 9 copies per gram of dry soil. Unexpectedly, fungal ITS copies did not significantly change along fire chronosequences across all soil layers regardless of fire severity (Table 1). Tukey's test also revealed no significant difference between fire stages. The relationships between fungal ITS copy number and pH, DOC, TC and NH 4 -N concentration were stronger in high-severity-fire sites than in low-severity-fire sites ( Figure S5).

Fungal community structure
Due to ITS1F and ITS4 being too long to pair-end, we analyzed them separately and received similar results ( Fig. 1 and Figure S2). Therefore, we used the ITS1F for further analysis. A total of 250,983 high quality sequences were obtained after denoising and quality filtering from ITS1F primer. We observed a total of 689 OTUs. The OTUs were classified into six phyla, including Ascomycota, Basidiomycota, Zygomycota, Chytridiomycota, Glomeromycota and Rozellomycota. Ascomycota was the Table 1 Results of linear mixed-effects models for the effects of fire stage on fungal diversity, fungal copy number and relative abundances of fungal guilds. Site is a random factor under low-and high-intensity fire. Significant P values are in bold The letter in parentheses indicated the best fitted model.  most dominant phylum, followed by Basidiomycota (Fig. 1). Archaeorhizomycetes was predominant across all samples with relative abundances ranging from 44 to 90%.
To better understand the differences in fungal responses to fire, we classified the fungi to different functional groups based on their life traits. A total 467 of 689 OTUs (approximately 16% of total sequences) were assigned to seven functional groups (Table 2). Among identified 467 OTUs, ~ 54% of fungal OTUs were saprotrophs, and ~ 29% OTUs were symbiotrophs. For the mycorrhizal OTUs, arbuscular mycorrhizal fungi were the most frequent and abundant, whereas Glomeraceae and Acaulosporaceae are the dominant families of all symbiotrophs. A total 242 of 689 OTUs (ca. 83% of total sequences) were assigned to unassigned groups due to its unknown ecological function. Among all unassigned groups, eighty OTUs (ca. 75% of total sequences) were identified as Archaeorhizomycetes.
Fungal community response to low-severity fire and other environmental factors Time since low-severity fire only had marginal effect on Shannon diversity in 0-5 cm layer, but not in other layers, as revealed by linear mixed-effect models (Table 1 and Fig. 2A). At 0-5 cm depth, Shannon diversity increased at 71d after lowseverity fire and then returned to pre-fire level by 2y ( Fig. 2A). At high-severity fire sites, Shannon diversity of 6-10 cm soil sharply decreased at 5y (Fig. 2B), and consistently declined over time in 11-20 cm soil ( Fig. 2C and D). Among all samples, the Shannon diversity had a weak positive relation to NH 4 − N (R 2 = 0.05, P = 0.001).
For low-severity fire, the PCoA plot showed that fungal community varied considerably along the time since fire gradient across all soil depths, and the variation seemed to be strongest in the surface layer (0-5 cm) (Fig. 3A-D). PERMANOVA, ANO-SIM and MRPP analysis showed that fungal communities significantly varied along the soil depths at each low-severity fire stage (Table 3, Figures S2A  and S2C). GDM analysis showed that time since fire was the strongest driver of total fungal community shifts with explaining 28.13% and 25.77% of the variance in 0-5 cm and 6-10 cm soil depths, respectively (Table 4), although there was a significant distance − decay relationship for total fungi in all low-severity fire samples ( Figure S6A). Below 10 cm soil depth, the geographical distance was most important driver that explained > 40% variance (Table 4). Besides time since fire and geographical distance, phenolics and water content played important roles in shaping fungal community (Table 4).
As time passed after low-severity fire, the relative abundance of Archaeorhizomycetes in 0-5 cm soil decreased significantly from 15 to 71d, and then increased at 2y, while saprotrophs showed an opposite trend ( Figures S7A and S7C). At 6-10 cm depth, Archaeorhizomycetes persistently declined but saprotrophs otherwise thrived ( Figures S7B and S7D). No significant patterns were observed below 10 cm soil. Except for a minor decrease at 15d, symbiotrophs were invariant to temporal variables across all soil depths according to linear mixed-effect model (Table 1).
Fungal community response to high-severity fire and other environmental factors Fungal communities appeared to progress over time after high-severity fire for every soil depth (Fig. 3E-H). PERMANOVA, ANOSIM and MRPP analysis showed that fungal communities did not vary along the soil depths at each high-severity fire  (Table 3, Figures S2B and S2D). The relative abundance of different fungal guilds did not show any significant pattern along the fire chronosequence (Table 1). Similar to low-severity fire samples, there was a significant geographical distance − decay relationship for total fungi among high-severity fire samples ( Figure S6B). Results of GDM analysis showed that total carbon and phenolics were the primary drivers, which explained 12.8% variance of total fungal community shifts in 0-5 cm, and 18.3% in 6-10 cm soil depths (Table 4). Shrub diversity explained less variation across all soil depths than soil properties and geographic distance (Table 4). Water content had no significant effect on fungal community across all sampling depths (Table 4).

Fire-responsive fungi
Within each fire severity, several frequent (≥ 50%) occurring OTUs significantly correspond to time lapses since fire ( Figs. 4 and 5). Classified as Mortierella genus, OUT6 declined in 6-15 cm soil over time after low-severity fire. Other OTUs (OTU9, OTU10, OTU11, OTU12, and OTU62), belonging to Archaeorhizomyces spp., exhibit diverse shift patterns along fire chronosequence. For example, OTU10 steadily declined with time since fire at depth of 6-15 cm. OTU62 increased across all soil depths. In the high-severity fire, OTU4, OTU11 and OTU12 belonging to Archaeorhizomyces significantly correlated with time since fire at 0-15 cm of soil, but not in the same way.

Responses to low-severity fire
Our results showed that the fungal copy number remained constant after short-term low-severity fire. ) represents the variance explained by fixed effects. Significant differences are reported as *, P < 0.05; **, P < 0.01; ***, P < 0.001. Fungal Shannon diversity for each fire stage at each soil layer is also given as the mean ± SE derived from 3 or 6 replicate measurements. Bars with the same letter were not significantly different in Tukey's HSD test reported from one-way ANOVA. Fire stage: 30 years post-fire (30y); 5 years post-fire (5y); 2 years post-fire (2y); 15 days post-fire (15d); 71 days post-fire (71d); undisturbed (UD) Some fungi that survive in spore banks (Glassman et al. 2016) might benefit from low-severity fire, promoting a rapid recovery. A decrease in soil pH after fire can favor some fungal growth (Pietikainen and Fritze 1995) and counteract the reduction of fungal copy number by heating. A linear relationship between the number of copies in the ITS region parsed by qPCR and the hyphal length of fungi can be established (Raidl et al. 2005). Therefore, to some extent, it is permissible to use fungal copy number to   Hence, the result also suggested that any changes in community structure along the fire chronosequence might be due to changes in individual relative abundances rather than to differences in the total number of fungal community among soils (Lee et al. 2017). Shifts in fungal Shannon diversity and fungal community composition along time since lowseverity fire gradient were pronounced in the surface (0-5 cm) layer (Figs. 2 and 3) coincides with other observations in Australia and the United States (Anderson et al. 2007;Brown et al. 2013;Sun et al. 2015). Although fungal Shannon diversity, relative abundance of Archaeorhizomycetes and saprotrophs in surface (0-5 cm) soil shifted after fire from 15 to 71d, they returned to pre-fire level at 2y ( Fig. 2A and Figure S7), indicating a high resilience to lowseverity fire. Low-severity fire decreased the relative abundance of Archaeorhizomycetes and increased saprotrophs in surface soil ( Figure S7). Previous studies suggest that Archaeorhizomycetes growth depends on root-derived materials (Rosling et al., 2011). Plant root derivatives are the main input of labile soil carbon (Fernandez and Kennedy 2016). Fire combusted the complete plant aboveground biomass in our study site (Flanagan et al. 2020), which may disrupt root derivative production and resulted in loss of Archaeorhizomycetes population. Low-severity fire along short duration increase the availability of fresh Fig. 4 Relative abundance of selected fire-responsive fungal OTUs along time since fire after low-severity fire. Fire stage: undisturbed (UD); 15 days post-fire (15d); 71 days post-fire (71d); 2 years post-fire (2y) organic matter and nutrients, which in turn stimulate the growth of saprotrophic fungi (Hartmann et al. 2014;Perez-Valera et al. 2019;Sun et al. 2015) in surface layers. Time since low-severity fire did not affect fungal composition and diversity in soil > 10 cm depth maybe because damp peaty soils are highly effective in delaying heat penetration (Davies et al. 2010), thus fungal communities found refuge below 10 cm depth during low-severity fire.
Although Archaeorhizomycetes dwindled after a burn episode, it remained the most abundant fungal class, accounting for 48%-74% of total sequences across all fire stages after a low-severity fire. Archaeorhizomycetes represents one of the most ubiquitous lineages of soil fungi and is an ancient class of fungi, yet only two species of this class have been identified (Menkis et al. 2014;Rosling et al. 2011;Rosling et al. 2013). Current knowledge showed that Archaeorhizomycetes might have K-strategies taxa with plant root-symbiosis (Hewitt et al. 2013;Naranjo-Ortiz and Gabaldon 2019) and markedly slow growth rates (Menkis et al. 2014;Rosling et al. 2011Rosling et al. , 2013. The persistence of dominant Archaeorhizomycetes may provide resilience in the ecological function of a system by alleviating post-disturbance lag times in ecological processes (Hewitt et al. 2013). Several OTUs (such as OTU9, OTU10, OTU11, OTU12, OTU62) belonging to Archaeorhizomycetes significantly varied as a function of time lapse since fire and with differential dominance patterns suggesting its non-linearity trait of responses. However, their ecological roles and interaction with environmental factors are unknown and need further studies.
A high proportional abundance of Archeorhizomycetes may overwhelm ericoid mycorrhizal fungi (ErMF) growth which may result in a low abundance of symbitrophs in this study (Table 1). Previous studies have shown that ericoid shrubs significantly attracted mainly Archeorhizomycetes (> 44%), which may suppress ErMF to a very minor contribution (< 0.1%) (Chronakova et al. 2019;Levy-Booth et al. 2019). Low abundance of ErMF (< 0.4% of total sequences) was detected at our sites and negatively related with Archeorhizomycetes (R 2 = 0.20, P < 0.05). As discussed earlier, the ecology of rootassociated Archaeorhizomycetes is largely unknown, so we can only speculate as to the role Archaeorhizomycetes play in those ecosystems and whether they actively interact with ericoid roots in assembling nutrients (Chronakova et al. 2019).
Responses to high-severity fire Meta-analyses indicated that Alaskan boreal forests would take at least a decade, or even up to two and a half decades, to completely shift back to the unburned state after a high-severity fire (Dooley and Treseder 2012;Holden et al. 2013). In boreal forests, most recent studies suggest that increasing burn severity induced the increase of fungal communities dissimilarity . In our sites, the highseverity fire killed all trees and led to 1-m peat loss. As a result, fungal community composition significantly differed from each other which suggests that at least several decades are needed for a complete recovery after high-severity fire in this evergreen shrub bog peatland.
In other studies following a high-severity wildfire, recovery of soil fungal community came after that of plant community because fungal recovery was more closely related to soil organic matter than the aboveground plant community (Dooley and Treseder 2012;Holden et al. 2013). In our study sites, dominant fire-tolerant shrub community recovered quickly (Table S1). Surprisingly, within the same high-severity fire site, depth have no effect on fungal community composition (Table 3) and fungal copy number (P = 0.13 and P = 0.77 for 5y and 30y, respectively), indicating high-severity fire eclipsed the depthdependent distribution patterns of the fungal community. It is important to note that 1 m of peat consuming high-severity fire may have selected for the well fire-adapted fungi (Reazin et al. 2016) across all studied soil depths. In addition, adaptation to soil nutrient environments is a niche process that can affect many microbial species due to the variance in their ability to break down recalcitrant organic matter (Talbot et al. 2013). Phenolics, one of recalcitrant organic matter, and total carbon increased with depth in the pre-fire site (R 2 = 0.23, P < 0.001; R 2 = 0.50, P < 0.001, respectively), but had no pattern in 5y and 30y samples. Long duration high-severity fire consuming substantial amounts of decay-resistance litter can release a pulse of nutrients (Jauhiainen et al. 2016;Sawyer et al. 2018), in turn, creating more similar soil chemical properties along soil profiles. Thus, within site scale, similar nutrient concentrations may form similar fungal community across soil depths.
Drivers of fungal community succession after low-severity fire Although geographic distance was an important predictor in fungal community assemblage (Talbot et al. 2014;Yang et al. 2020), low-severity fire events weaken spatial effects in fungal communities in surface layers (Table 4). Time since low-severity fire was identified in our study as the strongest driver for fungal composition turnover in the surface layer, which is consistent with reports in cold forest organic horizons (Yang et al. 2020). Time per se is not an environmental factor; rather, it integrates a range of changes, both biotic and abiotic, that occur over time (Cutler et al. 2017). In our studies, time since fire was closely related to soil pH, phenolics and water content in the top soil (Table S2). Progressive changes in these factors are likely to have impact on fungal community structure as well (Cutler et al. 2017).
Phenolics were shown to play important role in structuring fungal community in peatlands (Wang et al. 2021). Phenolics are potent inhibitors of 1 3 microbial breakdown, which regulate peat carbon sequestration via an'enzymic latch' (Freeman et al. 2001). Fungi are the major decomposer of polyphenol (Freeman et al. 2012), the community compositions of which being a key driver of an ecosystem's phenol-degradation capability (Waldrop et al. 2000). A decline in the abundance of inhibitory phenolics following fire may stimulate some fungi growth, which was suppressed by high concentration of phenolics before the fire (Freeman et al. 2012), in turn, reshaping the fungal community and potentially influencing carbon storage.
It now appears that pH is a unifying factor strongly shaping the community structure of bacteria and fungi (Fierer and Jackson 2006;Hartman et al. 2008;Glassman et al. 2017). At much broader spatial scales and wider pH ranges, pH was identified as the most important predictor of soil fungi (Tedersoo et al. 2014). Although pH showed a decline pattern since low-severity fire, pH was not identified as an important predictor of soil fungal community because pH variation in our study (3.2-4.6) are not as large as in other studies (4-> 8) (Fierer and Jackson 2006;Rousk et al. 2010;Griffiths et al. 2011;Tedersoo et al. 2014). Nevertheless, a decline of pH value favored fungal growth which could accelerate the recovery of fungi (Bárcenas-Moreno et al. 2011).
Water table level is of primary importance for the microbial community, as they are very sensitive to water and oxygen availability (Jaatinen et al. 2007;Kwon et al. 2013). In our undisturbed site, the mean water table (~ -17 cm) (Table S1) resulted in our sampling depths (0 − 20 cm) to be largely under aerobic conditions. The prescribed low-severity fire only slightly decreased the peat thickness and watertable level (Table S1), thus oxygen availability and the redox state of the surface peat were unlikely to be affected by water-table shifting (Etto et al. 2014), which in turn may have no significant impact on the fungal community after low-severity fire. Changes in water table position may however alter surface transpiration and evaporation, which contribute to the water content of the surface layers of peat (Waddington et al. 2015). In our study, water content of surface peat was decreased after low-severity fire, which had a significant impact on the fungal community (Table 4). Decreased water content could favor some complement of fungi which are able to degrade recalcitrant organic carbon to survive at low water regime environments. For example, the relative abundance of Dothideomycetes in the surface peat (0 − 5 cm), known to be involved in lignocellulose degradation , were significantly and negatively related to water table content (R 2 = 0.22, P < 0.005).
Drivers of fungal community succession after high-severity fire Not surprisingly, significant distance − decay relationship were also observed in high-severity fire samples, but substrate effects in our study confounded this relationship. Fungal community reassembly after high-severity fire was strongly linked to differences in post-fire soil properties which determined by sitelevel variation (Alcañiz et al. 2018;Certini 2005). Stratification patterns of the fungal community were diminished by high-severity fires. Soil properties (either phenolics or carbon) primarily shaped fungal community reassembly after high-severity fire across all soil depths. Furthermore, environmental heterogeneity may affect the distance − decay relationship for fungi (Green et al. 2004), which was the case in our study.

Limitations of the study
The main limitation is that the response patterns to both low-and high-severity fire were detected from different timelines. Missing multiple short-term time sampling events may have affected our detection of fungal community responses to different high-severity fires (Ludwig et al. 2018). For example, in highseverity fire, no considerable difference in fungal copy number between three temporal points (pre-fire, mid-term and long-term history) were detected, while the structure of fungal community varied among different times. Fungal copy number, a fungal biomass proxy, recover faster than fungal community structure, and thus the dynamics of fungal copy number may have been overshadowed due to our missing intermediate short-term sampling points. It has for example been shown that direct effects of microbes decimated by fire may be no longer viable one year after fire . High-severity burns may also result in an immediate reduction in fungal biomass (Holden et al. 2016) and shifts in community structure (Day et al. 2019). Unfortunately, few studies have concentrated on long-term (decades-centuries) changes in community composition and structure (Cutler et al. 2017;Yang et al. 2020), but they will be crucial for understanding the consequences of peat fire management under future climatic change conditions. Therefore, although complicated, studying long lasting post-fire fungal succession is important. Hence, we also analyzed fungal community response to high-severity fire in the longest-term fire site we could document (30y) in this study.
Another limitation was that we did not consider post-fire vegetation variation in living shrub biomass. In this study, both low-severity and high-severity fire had no significant effect on the relative abundance of symbiotrophic fungi (Table 1). For low-severity fire, this may be correlated with relatively stable shrub diversity (Table S1). Low-severity fires resulted in high tree survival and shallow depth of burn in which mycorrhizal colonization and diversity can persist (Pressler et al. 2019). For high-severity fire, the dynamics of symbiotrophs were not detectable due to missing intermediate time sampling points and shrub diversity recovered at 5y (Table S1).
In the light of the heterogeneity of soil properties and water-table levels between study sites even within same area, there is still high intra-variability (Alcañiz et al. 2018). More specifically researches need to pay attention to these spatial effects to learn more about their influence. In this way, it should be possible to identify the best prescribed fire practices for a specific area to manage and maintain carbon sequestration rates and community diversity.

Conclusion
Our study shows that time since fire, geographic distance and soil characteristic together explained a large proportion of variation in fungal communities, which is strongly dependent on soil horizons and fire severity. Fungal communities are highly resilient to low-intensity fires, while vulnerable to high-severity fire in our peatland. Hence, it is likely that the fungal communities recover within low-intensity fire interval in this shrub peatland. High-severity fire induced a significant community structure shift along the fire chronosequence and eliminated vertical distribution patterns over three decades. Pressler et al. (2019) also noted that had it not been for fungal resilience, future intensive fire events within three decades would further hinder the recovery of soil microbiome and the ecosystem process it regulates including plant communities recovery following fire. Taken together, these results suggest that soil fungi in some subtropical peatlands may have a high capacity to recover quickly following low-severity fires and, in turn, keep carbon stock stable. This study also improved our understanding on how low-severity versus high-severity fires impact soil fungal communities, important information, given the mounting size and frequency of fires in peatlands under climate change.