Introduction

Wildfire is a prevalent natural disturbance across the boreal forest region1, and its frequency and intensity are expected to increase with climate warming2. Forest fire directly affects soil carbon (C) cycling via immediate and high CO2 emissions from biomass combustion, and indirectly through long-term changes in ecosystem C dynamics and in soil organic matter chemical composition that occur during post-fire forest recovery and succession3.

Soil microbial communities play essential roles in forested ecosystems via nutrient cycling and decomposition of organic matter4. Bacteria are the most abundant group of soil microorganisms in terms of species5 and are among the first organisms to colonize dead wood after disturbance event and metabolize the easily accessible substrates6.

Fire dramatically reduces the biodiversity of the microbial community and biomass of soil microorganisms2,7,8. The effect of fire on soil microbial community is the greatest closest to the soil surface in the organic horizon and in the top few centimeters of mineral soil where heating is most intense9. Therefore, even low-intensity non-stand replacing fires may affect microbial communities through heat-induced mortality8,9. Soil microbial responses to fire are likely to change over the course of forest recovery after disturbance, and may only persist until aboveground plant communities regenerate10. Fire disturbances often cause changes in soil temperature due to loss of canopy cover and increases soil pH, both of which are likely to have large direct and indirect effects on soil microbial communities11,12. The increase in soil pH after fire favors bacterial growth13. The alteration of soil physical and chemical properties, like hydrophobicity, nutrient concentrations, and C quality by fire, may in turn have negative consequences for microbes7. Therefore, it is important to understand how microbial communities respond to the post fire environment. Fungi are considered to be important decomposers in forest soil and have been studied for fire-mediated changes14,15,16. Bacterial changes following fire have been studied more extensively than fungi12,17,18,19. However, investigations of bacterial dynamics over the course of a century and a half following fire are rare. There is also need to determine how fire affects bacterial communities and their ecosystem functions using molecular biological tools.

The advanced sequencing technology (e.g. Illumina MiSeq sequencing) allows us to better understand the environmental microbial community structure20. Functional gene arrays, such as the Geochip, are efficient for assessing the functional attributes of microbial communities in different environments21,22. The present study is part of a larger project investigating the microbial community structure and gene function across a 152-year boreal forest fire chronosequence. In this study, we characterized the responses of soil bacterial communities across three different periods of time post fire in a 152-year fire chronosequence (2-year, 60-year and 152-year after fire) by Illumina MiSeq sequencing, coupled with a high-throughput functional gene arrays, GeoChip 4.0. The aims of the study were 1) to evaluate the short- and long-term effect of forest fire and post-fire succession on soil bacterial communities; 2) to assess the functional potential of bacterial communities across the three times ranging from a few years after fire to over hundred years after fire. We hypothesized that the bacterial community after fire will differ with increasing time period since the last fire, and the community functions involved in different biochemical processes will change correspondingly.

Results

Information on the MiSeq sequencing

A total of 312 301 high quality sequences were generated across the nine soil samples after sequence de-noising and quality filtering, covering the three burned sites. The number of sequence per sample ranged from 27 699 to 42 617 with average of 34 700 ± 4 953 (mean ± SD) sequence. All the sequences were classified to domain bacteria.

Bacterial diversity after fire

A total of 3 692 OTUs (7 995 OTUs including singletons) were estimated across the three burned sites. The 2y site harbored the most OTUs of bacteria (2 369 OTUs) followed by 60y (2 027 OTUs) and 152y (1 728 OTUs) sites (Table 1). The 60y site showed the highest diversity (H′, 5.3) and evenness (0.72). The regression analysis indicated that the estimate of species richness (Chao 1 and observed OTUs) decreased over time since fire (P = 0.03, R2 = 0.74). No significant difference in diversity was observed among the three burned sites.

Table 1 Bacterial richness and diversity index for 16S region libraries in the soil from each site after fire.
Full size table

Bacterial community structure after fire

The three sites shared 32% of the 3 692 OTUs (1 192 OTUs) (Fig. 1a). The number of unique OTUs to each site ranged from 233 to 602, accounting for 37% of the total number of OTUs (Fig. 1a). The 2y site had the highest number of unique OTUs (602) and the 152y site had the least number of unique OTUs (233). The 2y and 60y sites shared most of the OTUs (48%, 1767 OTUs). The PCoA based on the relative OTUs abundance indicated that the three sites with different fire history formed distinct communities (Fig. 2a). The PCoA plot explained 43.1% of the observed variation, where the first axis explained 25.1% of the variations and separated the 152y site from the 2y and 60y sites. The second axis explained only 18.1% of the variation and separated the 2y and the 60y sites. Subsequent analysis of PERMANOVA confirmed that the community structure differed significantly among the three sites (P < 0.004). Interestingly, further analysis by separating abundant (>0.1% of total sequences) and rare (<0.1% of total sequences) OTUs indicated that the rare OTUs contributed to the difference in community structure between sites (Figure S1).

Figure 1
figure 1

Venn diagram showing the unique and shared OTUs from MiSeq sequencing (a), and presence/absence of functional genes identified by GeoChip (b) between the three burned sites. A total of 3 692 OTUs and 748 genes were detected across the three sites, respectively. Abbreviation: 2y: 2-year after fire; 60y: 60-year after fire; 152y: 152-year after fire.

Full size image
Figure 2
figure 2

Principle Coordinate Analysis (PCoA) plot (environmental variables as vectors) showing differences in bacterial community structure (a) and gene structure (b) with time since fire.

Full size image

The OTUs from all the samples were classified to 10 bacterial phyla (Fig. 3). The majority of the OTUs belonged to the Proteobacteria (41.3%), followed by Actinobacteria (15.4%), Acidobacteria (10.4%), Bacteroidetes (6.9%), Planctomycetes (5.8%), Armatimonadetes (3.9%), TM7 (2.2%), Verrucomicrobia (1.7%), Gemmatimonadetes (0.6%) and Firmicutes (0.2%) (Fig. 3). Protecobacteria was the most abundant phylum across all the samples (38.6%) followed by Acidobacteria (34.4%) and Actinobacteria (17.8) (Fig. 3). A small percentage of the sequences (4%) could not be classified to any phylum. The abundance of Proteobacteria increased over time since fire (Figure S2). The three most abundant phyla (Proteobacteria, Actinobacteria and Acidobacteria) accounted for more than 90% of the total sequence reads. The relative abundance of the detected phyla remained relatively stable in each site (Figure S2). No difference in the abundance was observed for each phylum between sites.

Figure 3
figure 3

Bar chart showing the phylum-level assignment for operational taxonomic units (OTUs) from three sites differing in fire histories as the relative proportions of OTUs and sequence reads.

Full size image

Among the 3 692 OTUs, 83% of them were assigned down to 33 classes, representing 95.4% of total sequence reads and 37% were assigned down to 132 genera, representing 52.3% of total sequence reads. Similar to phylum level, most of the classes did not show differences in the abundances between the sites. Only the class of Sphingobacteria and TM7 showed higher abundance at the 2y site than at the 152y and 60y sites (P = 0.03 and P = 0.04), respectively. Five genera showed significant differences in the abundance between the sites (Aciditerrimonas, Acidocella, Methylobacterium, Rudaea, Sphingomonas) and seven genera were represented only at the 2y site (Amnibacterium, Rhizobium, Spirosoma, Subtercola, Humicoccus, Kaistia, Marmoricola) (Table S1). Genus Sphingomonas had higher abundance at the 2y site than at the 152y site (P = 0.03). The abundance of Methylobacterium at the 2y site were also higher than that at the 60y site (P = 0.046, P = 0.03, respectively). The 152y site had higher abundance of Aciditerrimona and Acidocella than the 2y site (P = 0.034, P = 0.032, respectively).

Relationship between environmental factors and bacterial community structure

The soil bacterial community after fire was correlated with soil biogeochemical variables (R2 = 0.74, P = 0.01) (Table 2). The Distance Based Linear Model (DistLM) tests showed that the community structure was significantly correlated with soil temperature (P = 0.011), soil water content (P = 0.015) and soil pH (P = 0.017) (Table S2). Subsequent partitioning analysis showed that each of the three variables explained more than 20% of the observed total variation. Soil temperature and soil pH were positively correlated with the community at the 2y site, while soil water content was positively correlated with the community at the 152y site. The tree root biomass (TRootBiom), fungal ergosterol, soil C and soil N did not show significant correlation to bacterial community.

Table 2 The description of fire age characteristics and soil properties of the study sites a .
Full size table

Bacterial gene diversity and structure after fire

In total, we detected 62 860 probes that originated from bacteria in the GeoChip 4.0 in all samples, representing 748 functional genes (Table 3, Tables S3 and S4). These genes were classified into 11 functional categories (Table 3), including metal homeostasis (17 208), secondary metabolism (10 439), C cycling (9 443), sulfur (9 415), organic remediation (5 571), nitrogen (N) (2 598), stress (1 748), secondary metabolism (1 479), phosphorus (1 430), electron transfer (382) and others (3 153). Gene probes originating from fungi, archaea and virus were excluded from this analysis to match the bacterial taxonomic analysis. Among the 748 functional genes, 92% of them (688) were shared between the three sites (Fig. 1b). Only a few genes detected were unique to each site.

Table 3 Number of gene probes detected in each gene category from the sites with different fire histories.
Full size table

The highest gene diversity (H′) was observed at the 2y site followed by 152y and 60y site (Figure S3). The 2y site also had the highest number of detected probes (51 582 ± 432; mean ± standard deviation) (Table 3). The 60y site had the lowest gene diversity (H′) and the least number of probes detected (46 800 ± 2487; mean ± standard deviation). The gene diversity and the number of detected probes at the 2y site was significantly higher compared to the 60y site (P = 0.04, P = 0.03).

The PCoA analysis based on the intensity of detected probes revealed that the three sites formed separate gene profiles (Fig. 2b). Subsequent PERMANOVA analysis confirmed that the composition of functional genes between the sites differed significantly (P = 0.001 in all pairs). To visualize the pattern of detected functional probes among the sites, a cluster analysis was performed. The hierarchical cluster analysis revealed a difference in functional genes among the three sites (Fig. 4). The 2y and 152y sites were clustered together, indicating similar functional gene profiles, whereas both differed from the 60y-B site (P = 0.003, P = 0.003, respectively).

Figure 4: Hierarchical cluster analysis of functional gene probes in the three sites differing in fire histories based on the signal intensity of detected probes.
figure 4

The numbers after the site name represent replicates 1 to 3 at each site.

Full size image

Functional genes involved in different processes

Among the 11 gene categories, the number of detected probes in each category at the 2y site was significantly higher compared to the 60y site (Table 3). The signal intensity in each category was also significantly higher at the 2y site than that at the 60y site (Table S4). The genes involved in C, N and P cycling with significant difference in signal intensity between sites were shown in Fig. 5a–d.

Figure 5
figure 5

Normalized average signal intensity of genes involved in carbon cycling (a,b), nitrogen cycling (c) and phosphorus cycling (d) showing significant difference between sites. Data were presented as the mean ± standard deviation. The bars represent the standard deviations, and different letters in each panel represent Tukey’s significance at a P value of 0.05.

Full size image

The genes involved in carbon degradation (6 359 probes), carbon fixation (2961 probes) and methane oxidation and methanogenesis (123 probes) were detected in all samples. The genes involved in degradation of chitin (chitinase), lignin (glx, vdh), hemicellulose (ara), pectin (pme), starch (amyA, nplT) were detected at the 2y site with higher signal intensity than that at the 60y site (P = 0.05) (Fig. 5a). Higher signal intensity of several genes involved in carbon fixation of calvin cycle (FBPase, pgk and TIM), multiple systems (pcc) and reductive acetyl-CCoA pathway (codH, fthFs) were also detected with higher intensity at the 2y site compared to the 60y site (P = 0.05) (Fig. 5b).

For N cycling, genes involved in denitrification (1 135 probes), nitrification (471 probes), dissimilatory N reduction (389 probes), assimilatory N reduction (150 probes), N fixation (303 probes) and ammonification (110 probes) were detected in all samples and most of the genes showed similar signal intensity between the sites. Only the genes for N fixation (nifH) and denitrification (narG, nirS) showed higher signal intensity at the 2y site compared to the 60y site (P < 0. 05) (Fig. 5c).

In total, 1430 phosphorus cycling probes were detected, involved in polyphosphate degradation (1122 probes), polyphosphate synthesis (212), phytic acid hydrolysis (88) and phosphorus oxidation (8) were detected. However, only polyphosphate degradation gene (5f1_ppk2) and polyphosphate synthesis (ppk) gene showed higher intensity at the 2 y site than that at the 60y site (P < 0.05) (Fig. 5d).

Relationship between environmental factors and bacterial gene structure

Similarly to the taxonomic structure, the bacterial gene structure was also correlated to environmental variables (R2 = 0.81, P = 0.003) (Fig. 2b, Table 2). DistLM tests indicated the bacterial gene structure were significantly correlated to soil temperature (P = 0.031), soil water content (P = 0.022), soil pH (P = 0.026) and soil N (P = 0.005) (Table S2). Each of the variables explained more than 20% of the total observed variation. Soil pH and temperature were positively correlated with the gene structure at the 2y site, whereas soil water content and soil N were positively correlated to the gene structure at both 60y and 152y sites (Fig. 2b). The tree root biomass (TRootBiom) and soil C and did not show significant correlation to bacterial community.

Discussion

We investigated the bacterial community composition and gene structure across a 152-year chronosequence after forest fire. The recently burned site (2y) harbored the most bacterial OTUs with the highest number of detected gene probes. Both phylogenetic and functional structures of bacterial communities differed between the sites with different fire histories that were clearly correlated to the environmental factors. The pH and temperature affect mostly the community composition 2y after fire, whereas later soil moisture, soil N and N content and root biomass had greater impact to the bacterial communities.

The recently burned site (2y) site had the highest OTU richness and the richness decreased significantly over time after fire. Increased quantities of decomposable and burned material after the fire increase soil pH due to the production of K- and Na-oxidase, hydroxides and carbonates via ash deposition and accumulation23. The increased soil pH can favor bacterial growth over fungal growth13. Fire typically increases soil pH and other studies have shown decreases or no change in soil pH24. Similarly, the 2y site also harbored the highest number of gene probes, suggesting a positive correlation between biological richness and functional potential. Bacterial diversity has been reported to decrease with decreasing soil pH, in which a large variation of soil pH value was included25,26. However, we did not observe a difference in the diversity indices among the sites with different fire history. After the fire, a relatively small pH ranges (3.8–4.0) was evident, which might make it difficult to determine the impact on changes in diversity. In addition, the changing environmental factors after the fire led to a shift in community composition and resulted in a change in the evenness of each group. For example, due to a more even community composition, the diversity in the 60y old site was higher than that in the 2y old site, although the species richness was lower at the latter site. These could be possible reasons explaining no observed differences in diversity between sites.

As seen also in our study, it appears that environmental factors (e.g., soil temperature, pH, water content and soil N) have fundamental impacts on bacterial community structure7,27,28. A shift in temperature would lead to changes in microbial community structure and function29. Soil pH has been a key factor driving soil bacterial community across a large range of ecosystems25,30, and soil water content, organic C and N availability are also major determinants of soil microbial community composition31. Fire seems therefore to affect the soil bacterial community via altering the soil properties11,32. The loss of vegetation cover after fire disturbance can also contribute to other factors (i.e. loss of water, organic matter, etc.) and lead to higher soil temperatures in fire disturbed areas11,33. Higher soil temperatures also favor the enzymatic activity of microbes34. Our previous study conducted at the same sites showed that the fungal diversity was significantly affected by the fire and had long succession from 2 to 152 years after the fire16. These results indicate that bacterial community might recover more rapidly after the fire than fungi as fire does not seem to have any negative significant impact on bacterial diversity. Fire also reduced the biomass of trees and ground vegetation in our study sites and Köster et al.33 observed a significant reduction of fungal biomass 2y after fire in the same chronosequence. The decline in biomass was associated with an increase in the relative abundances of saprophytic fungi16. These findings support the conclusion that due to lower competition by mycorrhizal fungi, free living saprophytic organisms, both bacteria and fungi, become more abundant and dominate the communities after fire.

Many OTUs (1 381, ~37%) were unique to a specific site, suggesting an overall high beta-diversity. The sites with different fire histories formed distinct communities. Interestingly, the rare OTUs (<0.1%) significantly contributed to the community structural differences. The most abundant OTUs (>0.1%) showed similar distribution across the sites. This pattern has been observed also in other studies35 indicating that the dominant OTUs are not sensitive to environmental changes whereas the high number of OTUs with low abundance are critical when bacterial community responses are evaluated. Phylogenetically, Proteobacteria, Actinobacteria and Acidobacteria were the main phyla across the sampling sites, which accounted for 90% of the total sequences. Similar results were also observed in a Chinese boreal forest after wildfire36 and also in subalpine soil environments28. The relative abundances of bacterial phyla have earlier been shown to differ significantly from 4 and 16 weeks37 as well as from 1 and 11 years36 following a fire. In our study, however, there was no significant difference in the relative abundance of each detected phyla or class among the sites with different fire histories. The relative abundance of bacteria at high taxonomic scale remained its temporal stability. Only few genera showed differences in the relative abundance between sites, while most of them were present only in 2y site. These differences compared to previous studies could be due to the length of time between the fires as previous studies were conducted several weeks to 11 years after the fire occurred. In our study, we sampled much longer time after the fire occurred (2y to 152y).

The GeoChip data analysis demonstrated that the bacterial community gene structures were clearly separated among the three sites with different fire histories. The three replicates from each site clustered closely together indicating consistency of the biological replication. Overall, more than 90% of the detected genes probes were present in all samples. The minor remaining genes unique to each site suggested a low beta gene diversity in the study site. Notably, genes corresponding to nearly all major pathways of C (C degradation, C fixation and methane) and N (denitrification, nitrification, assimilatory and dissimilatory N reduction, N fixation and ammonification) were identified by GeoChip analysis, which might indicate large gene pool at all the sites irrespective the fire history. Microbial C cycling plays an important role in the forest ecosystem. The number of detected genes involved in C cycle was higher at the 2y site after the fire compared to older sites. Several genes (chitin, lignin, hemicellulose and pectin) involved in carbon degradation had higher signal intensity at the 2y site. Fire leads to the loss of readily degradable C and a corresponding decrease in soil moisture content38. The higher number and signal intensity of genes related to carbon degradation in recently burned sites might indicate that the potential for expression of carbohydrate-degrading enzymes was not affected or reduced despite the large reduction in C content caused by fire33.

Nitrogen is typically the limiting nutrient in northern forest soils39. Wildfire can cause long-term changes in N-cycling40. The nifH gene encoding nitrogenase reductase has been used as molecular marker for N fixation for a wide range of heterotrophic bacteria41. After the fire, the N content was reduced in forest soil, but the relative abundances for genes involved in N fixation (nifH) were not significantly impacted by fire12. On the other hand, Yeager et al.42 found the number of dominant nifH sequence types was greater in fire-impacted soils and Cobo-Díaz et al.43 also showed that the burned rhizosphere showed increased number of nifH gene copies. In our study, the nifH gene had higher signal intensity and number of probes detected in recently burned site compared to older sites. We propose that this finding indicate N fixation is involved in the compensation for N loss from the soil after burning.

Differences between young and older forests were observed in the gene probe number of the phosphorous cycling enzymes (Table 3). Fires are supposed to increase the availability of phosphorous due to mineralization7 and it seems that post-fire bacterial communities process the available phosphorous more actively. However, in litter bags analyzed at the same sites the phosphatases activities were higher in old sites44 which indicate that competitive replacement of bacteria by fungi affects also the function of bacterial communities. Altogether, the high activities of phosphorous cycling enzymes show that it is likely that early succession post-fire bacterial communities use actively the newly mineralized phosphorus.

The redundancies in the genes represented in the array compromised our ability to detect functional differences. The profiles of gene probes at all sites are diverse as the majority of them are present in all sites. Only a few differences are due to the presence/absence of certain genes or their abundance. It is important to bear in mind that the DNA-based GeoChip array reflects only the functional potential of microbial communities, but does not allow conclusions on how these genes are expressed and to what extend the actual function of the microbial population differs. To validate the functional process and population, other in-depth analysis (e.g., RNA-based transcriptome and functional activity assays) are needed. Similarly, most of the functional genes involved in organic C and N degradation pathways in our study were present in all sites. The activity of these genes and the bacteria involved need to be further validated.

In conclusion, the number of bacterial OTUs decreased over time since fire. However, fire did not show significant impact on bacterial diversity between 2- and 152-year after the fire. This observation differed clearly from the fungal diversity after fire at the same sites16. However, shifts in bacterial community structure and function among the sites with different fire histories were observed. Soil temperature, pH, water and C and N contents were the most important factors in shaping the bacterial community structures and function, suggesting that the fire-mediated changes in the soil biogeochemistry were strong drivers of the observed shifts in bacterial structure. Interestingly, the sites with different fire history formed separate bacterial gene clusters despite very long recovery time. This may indicate small differences in the potential to maintain essential biogeochemical processes in soil. This study provides functional insight on the impact of fire disturbance on soil bacterial community dynamics.

Materials and Methods

Study site and sample collection

The study site was located in a northern boreal subarctic coniferous forest of the Värriö Strict Nature Reserve in northeastern Finland (N 67°46, E 29°35). The detailed information on the sites and the soil sampling have been previously described16,33. Briefly, the site was dominated by Scots pine (Pinus sylvestris L.), with scattered downy birch (Betula pubescens Ehrh.) and Norway spruce (Picea abies (L.) Karst) trees. The ground vegetation consists mostly of bilberry (Vaccinium myrtillus L.), lingonberry (Vaccinium vitis-idaea (Lodd.) Hulten), black crowberry (Empetrum nigrum L.) and reindeer lichen (Cladina) species. The study was carried out at three sampling sites with different time since the last fire: 2 years (2y), 60 years (60y) and 152 years (152y) after fire (Table 2). The fire events on all the sites were not completely stand-replacing44.

At each sampling site, two plots were chosen with a distance of 100 m. Five consecutive soil samples with distance of 4 m each plot were collected, resulting in 10 samples per site. The humus layer samples (0.5 to 1.0 cm thick) were collected from a 0.25 m by 0.25 m quadrat of soil at after removing the litter33. The samples were mixed well and transferred to 1.5-ml Eppendorf vials. The samples were frozen at −180 °C in liquid nitrogen within a few hours after collection and transported to the laboratory in dry ice for subsequent DNA isolation. Continuous soil temperature and soil water content measurements were taken from the same location at 15-min intervals in the field. Soil temperature was measured using silicon temperature sensors (Philips KTY81-110, Philips semiconductors, Eindhoven, the Netherlands) and water content was measured using soil moisture sensors (Thetaprobe ML2x, Delta-T De-vices Ltd., Cambridge, UK) connected to a datalogger (DataTaker DT80, Thermo Fisher Scientific Australia Pty Ltd., Victoria, Australia). 10 soil cores (100 mm in length and 50 mm in diameter) from the same locations in each of the sites were also collected and transferred to the lab for the measurement of soil pH, soil C and N content, and root biomass in humus layer. The analysis of the soil properties in the laboratory have been previously described and published by Köster et al.33 (Table 2). Briefly, the soil humus layers were separated from the soil cores. Soil pH (H2O) of the humus horizons was determined using a glass electrode to stir the soil suspensions in demineralized water with a ratio of 10 ml of sample and 25 ml of water. The soil C and N content of the soil (including roots less than 2 mm in diameter) was measured with an elemental analyzer (varioMAX CN elemental analyzer, Elementar Analysensysteme GmbH, Germany) after the soils were dried in an oven at 105 °C for 24 h, sieved through a 2-mm sieve and ground with a ball mill (Retsch, Han, Germany). The soil C and N stocks were calculated for humus layers based on soil C and N concentration of the fraction and bulk density, respectively. For root biomass analyses, both the understory and tree roots were separated from the soil by washing, and sorted into living and dead fractions based on elasticity and toughness. The roots were identified as Scots pine, birch or other broad-leaved species, and understory (mainly dwarf shrubs and grasses) roots and rhizomes based on microscopic morphology and color. The ergosterol was extracted from soil samples by cyclohexane and 10% KOH in methanol. The amount of ergosterol was measured with high-performance liquid chromatography (HPLC) (HP Agilent 1100, Hewlett Packard, USA) using a Kinetex 2.6u C18 100A reverse-phase column 75 × 4.6 mm (Phenomenex ApS, Allerød, Denmark). All samples were collected in mid of August 2011.

DNA extraction, amplification of 16S rRNA gene, and Illumina MiSeq sequencing

Genomic DNA was extracted from 0.25 g (fresh weight) homogenized humus soil sample using the PowerSoil DNA isolation kit (Mo Bio Laboratories, Carlsbad, CA, USA), according to the manufacturer’s instructions. The genomic DNA was purified using the Gene Clean Turbo kit (MPBiomedicals, LLC, France), quantified with Qubit 2.0 fluorometer (LifeTechnologies, Eugene, OR, USA), and adjusted to a final concentration of 10 ng μl−1. Three individual extracted DNA in each site were randomly pooled together, representing three biological replicates per site (9 samples in total) and used for further Illumina MiSeq sequencing and GeoChip 4.0 microarray.

16S region V1-V3 was amplified in a 2-step PCR using the primers pA and pD45 containing partial TruSeq adapter sequences at their 5′ end, ATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT and GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT, respectively. The first PCR was done in two replicate 25 μl reactions using Phusion Hot Start II polymerase (Thermo Fischer) and cycling conditions consisted of an initial denaturation step at 98 °C for 30 s, followed by 15 cycles at 98 °C for 10 s, 65 °C for 30 s, 72 °C for 10 s, and a final extension for 5 minutes. After PCR the two replicates were combined and treated with Exonuclease I (Thermo Scientific) and Thermosensitive Alkaline Phosphatase (FastAP; Thermo Scientific). A second PCR was performed with full-length TruSeq P5 and Index containing P7 adapters and 1–5 μl from the first PCR as template. Cycling conditions were similar to the first amplification but with 18 cycles and 50 μl reactions with no replicates. Final purification was performed with Agencourt AMPure XP magnetic beads from Agencourt Bioscience (Beckman Coulter Inc, MA, USA). DNA concentration and quality were verified with Qubit (Invitrogen) and Bioanalyzer 2100 (Agilent), respectively. The final PCR fragments were pooled in equal concentrations and run on a MiSeq Sequencer (Illumina) using v2 600 cycle kit paired-end (325 bp + 285 bp).

Illumina MiSeq sequence data processing

The sequence data were analyzed using the mothur standard operation pipeline (SOP, Version 1.31.2)46. Briefly, pair-end reads were combined to contigs with minimum overlap 25 bp. Data were pre-processed to reduce sequencing and PCR errors, remove poor quality sequences, and preclustered with a distance of 2 bp using a pseudosingle-linkage algorithm implemented in mothur to minimize sequences that contain pyrosequencing errors47. All potentially chimeric sequences were identified using the mothur-embedded UCHIME program48 and removed. Unique sequences were aligned against SILVA database using the Needleman method49. The aligned distance matrices were clustered into operational taxonomic units (OTUs) using the average neighbor algorithm and 97% sequence similarity. OTUs with one sequence (singleton) were not included in downstream analysis. The sequences were classified to taxa using the RDP Naïve Bayesian rRNA Classifier tool version 2.0 with an 80% bootstrap confidence thresh-old50 on the RDPII and filtered to exclude chloroplast, mitochondria, eukaryotic and unknown sequences.

Bacterial richness was estimated by chao151 and diversity was assessed using Shannon indices (H′)52. To correct the difference in sample size, we used a randomly selected subset of 27 699 sequences per sample (minimum number of sequences recovered among all samples) to compare relative difference between samples. The OTUs were arbitrarily defined as ‘abundant’ when the relative abundances were above 0.1% and as ‘rare’ when the abundances below 0.1%. The analysis was done using mothur programme46.

All sequence data are available through the European Bioinformatics Institute (EBI) with project accession number PRJEB12433.

GeoChip hybridization and data processing

GeoChip hybridization and data processing were carried out by following the protocol described previously16. Briefly, 100 ng of genomic DNA from the triplicates in each of the three sites was amplified by rolling circle amplification using the TempliPhi kit (GE Healthcare) and a modified protocol22. Approximately 2 μg of amplified genomic DNA was mixed with random primers (3 μg/μl random hexamers; Life Technologies) and then labeled with 15 μl of labeling master mix (2.5 μl of dNTP [5 mM dATP, dGTP, and dCTP and 2.5 mM dTTP], 0.5 μl of Cy-3 dUTP [25 nM; GE Healthcare], 1 μl of Klenow fragment [40 U m1−1; Imer, San Diego, CA], 5 μl of Klenow buffer, 2.5 μl of water)53. The labeled genomic DNA was purified using the QIAquick purification kit (Qiagen), according to the manufacturer’s instructions, and then dried in a SpeedVac at 45 °C for 45 min (Thermo Savant). Labeled genomic DNA was hybridized on the GeoChip 4.0 microarray, as previously described21,53. The signal intensity of each spot was scored as positive and retained if the signal-to-noise ratio (SNR), calculated as (signal mean - background mean)/background standard deviation, was ≥2.0. Data normalization and pre-processing to remove poor quality spots (SNR < 2.0), normalization of the signal intensity of each spot based on the mean signal intensity across all genes on the arrays and transformation of the data using the natural logarithmic form were carried out as described previously21. The normalized signal intensities per gene were calculated as the sum of intensities of the probes per gene, divided by the total number of the probes detected in each gene, and averaged across the three replicates per sample54.

Statistical analysis

The bacterial OTUs were defined as the abundant OTUs (>0.1% relative abundance) and rare OTUs (<0.1% relative abundance). One-way ANOVA tests (Non-parametric Kruskall-Wallis) with Hodges-Lehmann estimate were used to identify differences in diversity and species richness between sites with different fire histories. The individual OTUs (Table S1) and gene categories (Table S4) that differ quantitatively between sites were also tested by ANOVA analysis. The bacterial community and gene structure was visualized as Principle Coordinate Analysis (PCoA) plot with Bray Curtis similarity using relative abundances of OTUs or normalized gene intensity in PRIMER v.6 34 with the add-on package of PERMANOVA+55. Prior to PCoA analysis the data were square root transformed to meet PCoA criteria. A PERMANOVA test was used to assess the significant difference in community structure between sites with different fire histories.

The process and functions used for the Geochip analysis were performed in the R56 gplots package57 for hierarchical cluster analysis (HCA). The gene diversity of each sample was estimated using Shannon’s index (H′).

To determine which environmental parameter (i.e. water content, soil pH, soil temperature, soil C and soil N) could significantly partition the variation in bacterial community and gene structure, a Distance Based Linear Model (DISTLM) analysis was used after PERMDISP procedure, and its significance was assessed by PERMANOVA (analysis of similarities; 999 Monte Carlo permutation tests). The values of the environmental parameters corresponding to each DNA samples (3 samples) in each site were calculated as average of the value of the three pooled subsamples, which was matching to each of the three DNA samples per site (Table 2).

Additional Information

How to cite this article: Sun, H. et al. Bacterial community structure and function shift across a northern boreal forest fire chronosequence. Sci. Rep. 6, 32411; doi: 10.1038/srep32411 (2016).