Soil Fungal Community in Norway Spruce Forests under Bark Beetle Attack

Bark beetle infestation is a widespread phenomenon in temperate forests, which are facing significant weather fluctuations accompanying climate change. Fungi play key roles in forest ecosystems as symbionts of ectomycorrhizal trees, decomposers, or parasites, but the effect of severe disturbances on their communities is largely unknown. The responses of soil fungal communities following bark beetle attack were determined using Illumina sequencing of soil samples from 10 microsites in a mature forest not attacked by bark beetle, a forest attacked by bark beetle, a forest destroyed by bark beetle, and a stand where all trees were removed after a windstorm. The proportion of ITS2 sequences assigned to mycorrhizal fungal species decreased with increased intensity of bark beetle attack (from 70 to 15%), whereas the proportion of saprotrophs increased (from 29 to 77%). Differences in the ectomycorrhizal (ECM) fungal community was further characterized by a decrease in the sequence proportion of Elaphomyces sp. and Russula sp. and an increase in Piloderma sp., Wilcoxina sp., and Thelephora terrestris. Interestingly, the species composition of the ECM fungal community in the forest one year after removing the windstorm-damaged trees was similar to that of the mature forest, despite the sequence proportion attributed to ECM fungi decreased.


Introduction
Bark beetle outbreaks represent the most important disturbance events in temperate and boreal forests [1] and, together with windstorms [2] and fire [3], cause severe damage to forest stands in the northern hemisphere.Extensive damage caused by mountain pine beetle (Dendroctonus ponderosae) and mountain spruce beetle (Dendroctonus rufipennis) were reported from Western Canada [4] and the United States [5][6][7].In European temperate forests, the European spruce bark beetle (Ips typographus) represents the dominant herbivore insect [3,8].
Bark beetle populations do not usually cause serious damage to the health of natural forests, but they can rapidly breed and spread as soon as conditions become favorable, e.g.sufficient breeding resources and habitats represented by fallen trees after a windstorm [3].These conditions have become increasingly relevant in a climate change context [1,7,8].Global warming and changes in precipitation patterns have resulted in the intensification of outbreaks [1,3] and can promote the survival, reproduction, and dispersal of insect herbivores [9], contributing to their higher pressure on forest ecosystems [1].
Soil fungi represent an important component of forest ecosystems.They are regarded as the dominant degraders of organic matter and mediators of carbon cycling [10].Ectomycorrhizal (ECM) fungi form a symbiosis with tree roots of several plant families including Pinaceae [11] and may influence plant fitness and plant community development [12,13].ECM fungi can also influence tolerance against herbivores [14] and promote recovery from herbivore infestation [15].
Even though fungi play such important roles in forest ecosystem functioning, information about their reaction to bark beetle infestation is sparse.The only study that has dealt with a bark beetle outbreak in spruce forest reported a decrease in mycorrhizal fungi while saprotrophic fungi increased due to the higher mass of dead plant material and gradual changes in the fungal community drive by substrate availability [16].Similar observations were made by Treu et al. [17] and Pec et al. [18] in pine forests destroyed by mountain pine beetle in Canada.
A massive European spruce bark beetle outbreak in the Tatra Mts.(Slovakia) spread several years after a severe windstorm in 2004 which destroyed 12,000 ha of spruce forests [19].The peak of subsequent bark beetle outbreak culminated in 2008 and 2009.The area of dead tree crown projections reached 90% in some forest parts [20], while other parts of these forests were less affected.Nevertheless, the first signs of infestation occurred in 2009 with tree mortality reaching almost 6% in 2012 [20].The Tatra Mts. were later affected by another windstorm event in 2014, which destroyed about 50 ha of spruce forest.
The aim of our work was to study how the trophic fungal groups were influenced by bark beetle infestation and identify key species that were most affected by this event.We compared the soil fungal communities from forests attacked by bark beetle to those from forests not attacked and with forests where trees were removed after a windstorm, which represents a commonly used management practice after first signs of bark beetle attack.We expected an overall increase in the proportion of saprotrophic fungi in forest damaged by bark beetle and a decrease in the proportion of ECM fungi accompanied by gradual changes in ECM fungal species composition in favor of stress-tolerant species in a forest attacked by bark beetle.A similar effect on trophic fungal groups proportions was expected in a forest with removed trees, but the increase in saprotrophs was expected to be not so distinctive, because of the lower levels of deadwood consisting mainly of branch piles.

Study Area and Sample Collection
The study area was located in the Tatra Mts., Slovakia (Figure 1), in originally unmanaged mature spruce forests with an admixture of larch (Lariceto-Piceetum), currently affected by bark beetle attack.All forest types were of similar age, ranging from 130 to 140 years, with canopy cover before disturbance events reaching from 45 to 60% and spruce proportion from 80 to 90% accompanied by larch.The altitudinal range is 1120-1240 m a.s.l., mean annual temperature 4.0-5.3• C, and mean annual precipitation 833-1000 mm; soil types are cambisols and podzols with predominantly granodiorite or moraine bedrock [20,21].
Four types of forest were selected: i) K1 = forest that was attacked by bark beetle, that had a 10% dead tree crown projection in 2015, and where bark beetle infestation started to spread in 2009; ii) K2 = forest that was destroyed by bark beetle in 2012 and where bark beetle infestation started to spread in 2007; iii) DES = forest with fallen trees removed after a windstorm in 2014; iv) REF = mature forest not attacked by bark beetle.K1 was a gradually decaying forest with still living spruce trees, admixture of larch and natural regeneration of spruce, fir, and rowan; Vaccinium myrtillus L. 1753, Oxalis acetosella L. and mosses prevailed in the understory.Deadwood was scattered throughout the microsites.K2 represented a forest completely destroyed by bark beetle in 2012 with standing or downed dead trees and poor natural regeneration, mainly spruce and rowan seedlings, with a dense herb cover dominated by Calamagrostis villosa (Chaix) J.F. Gmel., Vaccinium myrtillus, and Rubus idaeus L. The deadwood was more widespread around the microsites compared to K1. DES was a forest that was destroyed by a windstorm in 2014, where the fallen trees were removed one year before soil sampling, and where only piles of branches (i.e., logging debris) and a few standing spruce and larch trees remained.The regeneration was poor, dominated by spruce, larch, and rowan, with an understory dominated by Vaccinium myrtillus, Calamagrostis villosa, and Avenella flexuosa (L.) Drejer.The mature forest (REF) consisted of spruce with larch admixture, spruce, and larch regeneration and an understory dominated by Vaccinium myrtillus, Avenella flexuosa, and Calamagrostis villosa.The deadwood was sparse and mostly in the form of branch piles.Soil samples were collected in September 2015.In each forest type, 10 microsites were selected from an area of 100 × 100 m along with two approximately 140 m long crossed diagonals that bisected the center of the sample area.The distance between microsites (approximately 30 m) was inferred from the observed variability of particular patches in mountain spruce forest.Five soil cores (3.5 cm in diameter, 20 cm depth) were collected per each microsite in a 5 m diameter circle and dominant living trees, herbs, and type of dead wood were recorded.The tree, herb, and dead wood cover was estimated in a 5 m diameter circle using a visual assessment based on crown/herb leaf/dead wood surface projection.The dead wood was categorized into four classes (lying logs, hanging logs, piles of branches, and uprooting) and had a threshold diameter between branch and log of 7 cm.The soil cores were taken from places without piles of debris or litter accumulation.The litter-and mineral soil horizons were then removed, and the material from organic soil horizon was passed through a 5 mm sterile sieve.The fresh soil from organic horizon from five cores per each microsite was then mixed together and stored at -20 • C.

Soil Analysis
Approximately 30 g of organic soil material were used for soil analysis.C and N concentrations were measured using a Flash 2000 CHN analyzer (Thermo Fisher Scientific, Waltham, MA, USA).The pH (H 2 O) of the soil suspension was determined by the potentiometric method on an OP-208 pH-meter (Radelkis, Budapest, Hungary) equipped with a combination glass electrode.

DNA Extraction, Amplification, and Sequencing
Total DNA was extracted from ~300 mg of soil using the NucleoSpin Soil Kit (Machery-Nagel Gmbh & Co., Düren, Germany) in three replicates per microsite.The extracted DNA from three replicates was then pooled and used as a template for ITS2 amplification (gITS7/ITS4 primer combination; [22,23]).The PCR reaction was performed in three replicates in a 25 µL reaction mixture, containing 1× reaction buffer (Bioline, London, UK), 1.5 µl BSA (10mg/ml), 1 µL of each primer (0.01mM), 0.75 µL of a polymerase mix (4U:96U of Pfu DNA polymerase (Fermentas, Waltham, MA, USA) and MyTaq polymerase (Bioline, USA)), and 1 µL of template DNA filled with dH 2 O up to 25 µL.The amplification conditions were 94 • C for 5 min, followed by 35 cycles of 94 • C for 30 s, 56 • C for 30 s, 72 • C for 30 s, and a final extension of 72 • C for 7 min.Three PCR replicates were pooled and purified using the MiniElute PCR purification Kit (Qiagen, Hilden, Germany).The concentration of purified PCR amplicons was quantified using a QUBIT 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA).The sequencing of fungal amplicons was performed on a MiSeq (Illumina, San Diego, CA, USA) at the Institute of Microbiology of the CAS, v. v. i.

Data Processing and Analysis
The sequence data were processed using pipeline SEED 1.2.3 [24].Chimeras were removed using USEARCH 7.0.1090[25] and the ITS2 region was extracted using ITS EXTRACTOR 1.0.8[26].Sequences were then clustered using UPARSE at a 97% similarity threshold, and consensus sequences were created using MAFFT 7.215 [27].The consensus sequences were then identified using a reference database for fungal ITS classification based on UNITE [28] release 7.1 2016-11-20 (plus plant and animal sequences from GenBank).To maintain a comparable sample size, 2800 reads per sample were randomly selected.Because the operational taxonomic units (OTUs) represented by a small number of reads were found to be not reproducible [29], only OTUs represented by more than 1% of read counts per sample were included in further analysis.After discarding non-reproducible and non-fungal clusters, 184,840 reads remained, representing 229 OTUs.The assignment of OTUs at the species level was based on the 97% similarity threshold, to the genus level on at least an 85% similarity threshold, for a coverage higher than 70%.The Shannon-Wiener index, species richness, and evenness were calculated using the pipeline SEED.The assignment of fungal species to trophic groups was done according to [30][31][32][33][34][35][36][37][38][39][40][41][42][43].Sequence data were deposited in the GenBank database (dataset number SRP157877).

Statistical Analysis
Analysis of variance in program R 3.4.3[44] was used to test the differences in values of organic soil pH/H 2 O, carbon concentration (C), nitrogen concentration (N), and the C/N ratio among the four forest types.Similarly, differences in fungal trophic groups (ECM fungi, saprotrophic fungi, parasitic fungi) among the forest types were tested.Logarithmic transformation was used in case of parasitic fungi data.
Variation in the fungal communities of the microsites, based on the number of sequences, was analyzed by multivariate methods in Canoco 5 [45].We analyzed three datasets: fungal sequence abundances, ECM fungal sequence abundances, and saprotrophic fungal sequence abundances, all logarithmically transformed.A principal component analysis (PCA) was first conducted to determine the total variability in the datasets.Redundancy analyses (RDAs) with forward selection of environmental variables with significant effect (Monte Carlo permutation test, 999 permutations) were then conducted for each dataset.As environmental variables were tested, coverage percentages of the tree layer, dominant tree species (Picea abies (L.) H. Karst., Larix decidua Mill., Sorbus aucuparia L.), dominant species of herb layer, including dwarf shrubs (Calamagrostis villosa, Vaccinium myrtillus, Vaccinium vitis-idaea L., Rubus idaeus, Avenella flexuosa, Epilobium angustifolium L., Oxalis acetosella, Dryopteris sp., Luzula sp., and mosses), litter, and dead wood.The forest type and soil characteristics (pH/H 2 O, C(%), N(%), and C/N) were also tested.The effect of forest type was also tested in separate analyses.Environmental variables were also analyzed in a PCA standardized by species, which shows how the forests differ in their environmental characteristics.

Fungal Community
In total, we detected 229 OTUs: 92 were mycorrhizal species, 79 were saprotrophic species, 12 were parasitic species, and 3 were endophytic species, and for 43 of them, either they were not identified at the genus level or available information about trophic status was insufficient.The diversity and evenness were similar in all forest types, with a Shannon-Wiener diversity index of ~3.5 and an evenness of ~7.3.Species richness was slightly higher in the K2 site compared to K1, DES, and REF (Table 1).Mycorrhizal relative sequence abundance inferred from OTU sequences belonging to mycorrhizal species was significantly lower (p < 0.05) in K2 compared to K1, DES, and REF (Table 1, Figure 2a, and Table S1b).A significantly lower number of mycorrhizal sequences was detected also in DES compared to REF (p = 0.05).The opposite trend was found in the case of saprotrophic fungi with a significantly higher (p < 0.01) relative sequence abundance in K2 compared to K1, DES, and REF (Table 1, Figure 2b, and Table S1b).Parasitic fungi were most abundant in K2 but non-significantly, while REF contained a significantly lower (p < 0.05) amount of parasitic fungi sequences compared to the other sites (Table 1, Figure 2c, Table S1b).The total number of OTUs assigned to ECM fungi, as well as saprotrophic and parasitic fungi, is listed in Table 1, including dominant ECM and saprotrophic species.The differences in fungal communities among forest types (Figure 3) were caused mainly by ECM orders (Figure 4), with Atheliales prevailing in REF (16% of sequences on average), but being substantially reduced in K2 (3%), along with Russulales (K1 18%, K2 1%, DES 11%, REF 19%).In contrast, Mortierellales containing saprotrophic species became more abundant in K2 (22%) and DES (13%), but less abundant in K1 and REF (6%).Eurotiales, including ectomycorrhizal, saprotrophic, and parasitic species, prevailed in DES (11%) and REF (12%) but decreased in K1 (6%) and K2 (2%).The PCA ordination diagram showed that the forest types were well separated based on fungal species composition (Figure 5a).Axis 1 accounted for 12.9% of the variability and followed the difference between mature forest (REF) and forest destroyed by bark beetle (K2), with DES overlapping both.The bark beetle-attacked forest (K1) was separated on Axis 2 from the other three forest types (10.4% of explained variability).The forest types were well separated according to ECM species composition with tiny overlaps between K1 and REF, DES and REF, and K2 and DES (Figure 5b).Similarly, to fungal composition, Axis 1 largely depicted the difference between K2 and REF, Axis 2 separating out K1.The forest types were also well separated based on saprotrophs, with only a small overlap between REF and DES (Figure 5c), and similar pattern as previous diagram depicting fungal composition.Microsites in mature forest (REF) are more variable in terms of ECM species composition compared to the saprotrophic ones, whereas microsites in forest destroyed by bark beetle (K2) show an opposite trend (Figure 5b,c).

Vegetation
Tree cover was the highest in K1 and REF forests, being significantly greater than in the more heavily disturbed K2 and DES (Table 2).The understory was dominated by Vaccinium myrtillus in all forest types, accompanied by Oxalis acetosella in K1, Callamagrostis villosa in K2 and DES, and Avenella flexuosa in REF (Table 2).Spruce seedling regeneration was the highest in REF, whereas almost no regeneration was found in the other forest types (Table 2).A high cover of dead wood mass was found in DES in the form of piles of branches and K2 in the form of lying logs.These were present to a lesser extent in K1 and REF, but piles of branches were also present at several REF microsites (Table 2).For detailed information, see Table S1a.

Soil Properties
The organic soil was strongly acidic in all forest types, with a mean pH of 3.9 ± 0.3 (Figure 6a).Carbon and nitrogen concentrations were the highest (p < 0.05) in K1, whereas the lowest concentrations were in DES (Figure 6b,c).The highest C/N ratio was detected in K1 and DES with means of 25 ± 0.7 SE and 25 ± 0.5 SE, respectively (Figure 6d).The K1 forest had the highest within-microsite variability in all measured soil properties.

Vegetation and Soil Properties vs. Fungal Communities
Based on the redundancy analysis (RDA) with explanatory environmental variables, forest type explained 18.7% of detected variability within the fungal community composition.Excluding forest type as an explanatory variable, tree cover, the coverage of Vaccinium myrtillus, Oxalis acetosella, and Drypteris sp., the percentages of soil nitrogen and carbon, and latitude (N • ) explained 22.2% of the variability among the microsites (Table 3).The ECM community was mainly influenced by forest type, accounting for 13.4% of the explained variation and further by tree cover, cover of Dryopteris sp., and latitude (N • ) (Table 3).Saprotrophic species were affected by forest type accounting for 14.0% of the explained variation.Excluding forest type, tree cover, the percentages of soil nitrogen and carbon, pH/H 2 O, and the presence of Sorbus aucuparia and Calamagrostis villosa explained 11.8% of the variability (Table 3).Correlation of the environmental variables with each forest type is shown in Figure S1.The most important variables are tree cover, carbon and nitrogen concentrations, and cover of Rubus ideaus and Oxalis acetosella.

Discussion
Our data confirmed the expected hypothesis that tree dieback will cause changes in proportions of sequences and community composition of soil fungi.These changes were primarily driven by shifts in the occurrence of ECM fungi, because they are directly dependent on carbon provided by their host trees [14].
The decrease in ECM relative sequence abundances was apparent within one year after tree removal (DES), but surprisingly species composition remained comparable with the mature forest (Figure 4).Though it is not possible to determine in absolute terms how the abundance of the ECM species was affected, the fungal community composition in forest with removed trees (DES) was comparable to that in mature forest.This indicates the ability of many ECM species to survive for some time in the absence of living trees.The possibility that we had DNA from dying mycorrhizae or spores is quite low, because dying mycorrhizae would not persist long [46].Similarly, spores in the soil would have a different composition than the living fungal community, because the detected fungi in the mature forest differ in their investment to fruit bodies and hence expected spore production [47].
The direction of ECM species change is in agreement with our data from stands affected by windstorms [48].Microsites in the mature forest (REF) and the forest with removed trees (DES) were both dominated by Elaphomyces sp., different Russula species, and Tylospora asterophora, whereas the K2 forest, destroyed by bark beetle, was dominated by Piloderma species, Inocybe soluta, Wilcoxina sp., Thelephora terrestris, and T. asterophora; T. terrestris, Piloderma sp., Tylospora fibrillosa, T. asterophora, and Wilcoxina sp. also dominated in stands 10 years after fire and windstorms [48].Buscardo et al. [49] and Ford et al. [50] considered T. terrestris and Wilcoxina species to be stress-tolerant, being able to survive disturbances or rapidly re-establish after disturbances [51].Propagules of T. terrestris and Wilcoxina are also commonly found in spore banks [52].The ecologies of Piloderma and both Tylospora species are not clear, which could be attributed to taxonomical difficulties and the presence of more strains differentially adapted to local conditions [53].Tedersoo et al. [54] suggested that they contribute to wood degradation, enabling them to compete with saprotrophs for resources.The ecology of the Russula genus is also uncertain.In our case, it dominates in the healthy mature forest (REF) and forest with removed trees (DES), whereas Štursová et al. (2014) [16] treated it as potentially saprotrophic being able to survive bark beetle outbreak.This difference is most likely caused by different dominant species in forests in the Tatra (R. mustelina, R. puelaris, R. vesca) and Sumava Mts (R. ochroleuca Pers.1796), which emphasizes the need to take into account species identity in any ecological interpretations.The bark beetle-attacked K1 represented a transitional forest stand with species belonging to both the mature forest (e.g.Russula species, Elaphomyces sp.) and destroyed stand (e.g.T. terrestris), which resulted in a temporal increase in species diversity (Table 1).Relative abundances of ECM fungi in K1 were similar to that of REF forest, supporting the hypothesis that an ECM fungal network may be stable if enough trees remain alive [55,56].It is not yet clear if leaving such attacked forests for natural development would help maintain a high species pool in the landscape compared to forests being immediately cut after a first sign of bark beetle occurrence.Generally, ECM fungal communities undergo transitions after disturbances, shifting from species with higher demands on resource availability to species more adapted to unfavorable site conditions.ECM fungi compete for resources and, as the resource availability becomes diminished, the weaker competitors are excluded from the stand [57].
In comparison with changes in the ECM communities, the qualitative change in saprotrophic communities was not as straightforward.We expect that this is due to their lower dependence on living trees and smaller body size, enabling them to easily find their optimal niche, even though they were found to be more susceptible to soil abiotic factors compared to ECM fungi (Table 3).Rhodotorula sp. was the dominant genus in all forest types.Together with Cryptococcus sp., which was slightly decreased in K2, it represents stress-tolerant yeast organisms [58].Cryptococcus sp. has also been found to increase its abundance after bark beetle infestation [16].Mycena spp.increased in both bark beetle-affected forest types, with M. pura increasing in K1 and M. galopus in K2.Mycena species are free-living filamentous fungi, which compete with yeasts for resources [58].Their increase after bark beetle infestation could be attributed to higher resource availability represented by dead wood and decaying root systems.Rhizidium phycophyllum decreased only in the mature forest.Even though it is considered to be a saprotroph, it is also described in association with green algae [59].The reason why it is abundant in destroyed forests is not clear.Due to our focus on soil fungi, we were only able to detect wood-decaying saprotrophs by accident, because it takes decades for attacked trees to fall following a bark beetle outbreak, so the delayed availability of the substrate will contribute to extended delays in changes in saprotrophic composition [18].

Conclusions
Bark beetle outbreak significantly influences soil fungal communities.ECM fungal communities in forests affected by bark beetles are dominated by stress-tolerant species and species forming fruit bodies preferentially on wood.It remains unclear how long species of mature forest persist in stands affected by such events.Similarly, there is the question of whether unmanaged forests attacked by bark beetle, harboring both species of mature forest and stress-tolerant species, could contribute to future biodiversity and regeneration as a potential reservoir of ECM propagules.We expect that the weaker differentiation of soil saprotrophs was due to their lower dependence on living trees and smaller body size, making it easier for them to find their optimal niche.

Figure 3 .
Figure 3. Average relative abundance of sequences assigned to fungal orders from 10 microsites per forest that is attacked by bark beetle (K1), destroyed by bark beetle (K2), with removed trees (DES), and mature (REF).

Figure 4 .
Figure 4. ECM fungal taxonomic composition assigned to ECM orders as noted by average absolute sequence abundance from 10 microsites per forest that is attacked by bark beetle (K1), destroyed by bark beetle (K2), with removed trees (DES), and mature (REF).

Figure 5 .
Figure 5. Principal component analysis (PCA) summarizing variation in the fungal community based on the number of sequences from 10 microsites per forest that is attacked by bark beetle (K1), destroyed by bark beetle (K2), with removed trees (DES), and mature (REF): (a) fungal sequence abundances; (b) ECM fungal sequence abundances; c) saprotroph sequence abundances.

Table 1 .
Species richness, the relative sequence proportion, the number of species and dominating ECM, and saprotrophic species in the forest.K1 = attacked by bark beetle; K2 = destroyed by bark beetle; DES = with removed trees; REF = mature and not attacked by bark beetle.
*average relative sequence proportion in % and standard error from 10 microsites per forest type.

Table 2 .
Mean cover of trees, dominant herbs, spruce regeneration, and dead wood mass in forests: K1 = attacked by bark beetle, K2 = destroyed by bark beetle, DES = with removed trees; REF = mature.
*average cover in % and standard error from 10 microsites per forest type.

Table 3 .
Variation in the fungal community based on the number of sequences using RDA; environmental variables used as explanatory variables: a) forest type as an explanatory variable, b) environmental variables excluding forest type, and c) all environmental variables (N%, nitrogen content; C%, carbon content; N • , north latitude) in forests: K1 = attacked by bark beetle; K2 = destroyed by bark beetle; DES = with removed trees; REF = mature.