Reindeer grazing alter soil fungal community structure and litter decomposition related enzyme activities in boreal coniferous forests in Finnish Lapland

Reindeer grazing in northern boreal zone affects forest floor vegetation heavily and alters the vegetation structure. However, the effect of grazing on soil fungal communities, which are intimately linked to plants, is not currently known. Therefore, our objectives were to investigate changes caused by reindeer grazing on soil fungal communities, litter decomposition rate and litter degrading extracellular enzyme activities. The study was conducted in four areas divided into grazed and non-grazed sites (all together 38 sample plots) in northern boreal forests in Finnish Lapland. Fungal communities were analyzed from humus with high-throughput sequencing technology (454-pyrosequencing), and litter mass loss and extracellular enzyme activities were analyzed after a one-year litterbag experiment. The results showed that grazing significantly affected the fungal community structure and the abundance of certain fungal genera and species. Grazing also decreased laccase and enhanced cellobiohydrolase I activities from the litterbags. Our study is one of the first to describe detailed fungal community composition in sites with long-term history of reindeer grazing and exclusion. Our results indicate that reindeer grazing alter fungal community structure and litter degradation related enzyme activities in the northern boreal forest soils.


Introduction
Semi-domesticated reindeer (Rangifer tarandus L.) are the most numerous large mammalian herbivores in the northern boreal zone, especially in areas where their husbandry is practiced. Currently, there are on average 200 000 reindeer (number in winter after culling) (Suominen and Olofsson 2000) grazing freely in Finnish Lapland, on an area covering one third of Finland. The reindeer density is on average 1.5 individuals per km 2 (Susiluoto et al., 2008;Turi, 2002), although the density varies between reindeer herding districts and are above the ecosystem's carrying capacity in many areas in Finnish Lapland (Suominen andOlofsson 2000, Kumpula et al. 2014).
Reindeer are strongly affecting the forest floor vegetation in northern latitudes and their direct and indirect consequences on ecosystems are important to understand. Compared to many other large mammalian grazers, reindeer feed on only a minor part of net primary production (NPP) but alter the vegetation structure by selective grazing and trampling. The main impact of grazing on forest floor vegetation is the reduction of lichen biomass (den Herder et al., 2003;Köster et al., 2015b;Stark et al., 2000;Susiluoto et al., 2008). A reduction of tree regeneration, total ground vegetation biomass, mosses and vascular plants in forested ecosystems has also been observed (Köster et al., 2013(Köster et al., , 2015bOlofsson et al., 2010). Furthermore, reindeer affect plant diversity and coverage, change soil microclimate, moisture conditions, nutrient cycling and organic matter decomposition (Akujärvi et al., 2014;Köster et al., 2013Köster et al., , 2017Köster et al., , 2015bOlofsson et al., 2004;Stark et al., 2010;Vowles et al., 2017;Väre et al., 1996).
The less is known about the potential effects of reindeer grazing in the boreal forest ecosystem. Grazing reduced lichen mat cover resulting in faster soil warming in spring and higher soil temperatures during growing season, which can enhance Scots pine (Pinus sylvestris) growth (Fauria et al., 2008). On the other hand, reduced lichen mat can increase soil moisture fluctuation and decrease plant growth (Olofsson et al., 2010). Grazing can also decrease soil microbial respiration and activity in boreal forests humus, which was suggested to result in decreased decomposition rate (Väre et al. 1996, Stark et al. 2003. In Alaskan taiga, mammalian browsers (moose and snowshoe hares) were observed to significantly reduce the annual production, decomposition and survival of fine roots (Ruess et al., 1998), and to reduce ECM infections of fine roots of willow and balsam poplar compared to exclosures protected from browsing (Rossow et al., 1997). However, in Finnish northern boreal forest, reindeer grazing was not observed to affect tree root biomass (Köster et al., 2013(Köster et al., , 2015b and field layer respiration was found to be higher on grazed sites compared to nongrazed sites due to higher soil temperatures (Köster et al., 2017, Köster et al., 2018. Reindeer lichens (Cladonia spp.) have been reported to produce allelopathic extracts that can inhibit the growth of many mycorrhizal fungi, including mycorrhiza of Scots pine and dwarf shrubs in pure and axenic synthesis cultures (Brown and Mikola, 1974). Since mycorrhiza formation is essential for nutrient uptake of woody plants in boreal ecosystems (Smith and Read, 2008), growth of pine seedlings may be positively affected by reindeer grazing due to reduction of lichen cover. Lichen extracts were found to reduce seed germination and mycorrhizal colonization percentages of roots on pine seedlings in a greenhouse experiment (Sedia and Ehrenfeld, 2003). In contrast, Kytöviita and Stark (2009) demonstrated that in microcosms pure usnic acid, which is the most abundant secondary metabolite of the lichen Cladonia stellaris, or lichen fragments did not retard the growth of pine seedlings or negatively affect the plant nutrient uptake. In addition, usnic and perlatolic acids produced by C. stellaris lichens showed no evidence of antimicrobial activities, as these substances did not leach out from the lichen thallus to the rainwater, inhibit microbial respiration or enrich in soil underneath the lichen mat (Stark et al., 2007). However, to our knowledge, no field studies under natural conditions have been carried out to confirm or contest these findings in lichen rich boreal forests.
Only a few published studies have investigated the effects of reindeer grazing on soil microbial communities in northern boreal ecosystems. Generally, these studies have used phospholipid fatty-acid analysis (PLFA) and concluded that grazing had no clear effect on microbial community structure or biomass (Bardgett et al., 1996;Francini et al., 2014;Rinnan et al., 2009;Stark et al., 2010Stark et al., , 2008. As PLFA is not a highly sensitive method to describe microbial, and especially fungal community changes, current methods based on new sequencing technologies may therefore better reveal these patterns. Vowles et al. (2018) studied mycorrhizal communities in Scandies Mountain Range with ingrowth bags and found that although grazing did not affect the diversity of mycorrhizal species, it decreased extramatrical mycelial biomass and the abundance of genus Cortinarius in mountain birch forest site. Thorough, sequencing-based analysis on the effects of reindeer grazing on soil fungal communities are still lacking, especially from the northern boreal forests.
With variety of extracellular enzymes, fungi are the predominant decomposers of plant residues and soil organic matter (SOM) (Lundell et al., 2014;Rytioja et al., 2014). In mesotrophic tundra heat in northernmost Norway, microbial respiration, ß-glucosidase, acidphosphatase and leucine-aminopeptidase activities were higher and Nacetylglucosaminidase activity was lower in heavily grazed areas compared to the lightly grazed areas (Stark and Väisänen, 2014). Sedia and Ehrenfeld (2006) tested litter decomposition rates and soil enzyme activities under lichen and moss covers and found that the decomposition rate tended to be slower underneath lichen than moss cover.
In the present study, our objectives were to describe the fungal community structure and to estimate the overall, long-term impact of reindeer grazing on soil fungal communities in grazed and non-grazed sites in northern boreal forest. Furthermore, our objectives were to show how grazing affects microbial communities' functionality by comparing the litter decomposition rates and extracellular enzyme activities between the grazed and non-grazed sites. To study these aims, we had four unique natural forest areas separated by 20, 55 or 100 years old fences, creating sites with and without reindeer grazing side by side. This allowed us to test whether in long-term reindeer grazing as a whole is a determining factor that shape the soil fungal communities in natural environment. We hypothesized that vegetation changes caused by reindeer grazing alter fungal community structure and accelerates the decomposition rate of litter. We also hypothesized that ECM fungi would be more abundant in the grazed than non-grazed sites, due to the allelopathic substances secreted by Cladonia lichens (Brown and Mikola, 1974). To our knowledge, this is the first study to document long-term changes in fungal communities caused by reindeer grazing in northern boreal forest by using high-throughput sequencing.

Sites and sampling
Our study areas were located in northern boreal coniferous forests at close vicinity of ICOS/SMEAR I station in Värriö (67°46′N, 29°35′E) and the Arctic Research Centre of the Finnish Meteorological Institute in Sodankylä (67°21′N, 26°38′E). We had four sampling areas: Nuortti 1, Nuortti 2 and Kotovaara in Värriö and one area in Sodankylä. All areas were non-managed, nature reserve forests with the average age of the dominant trees in Nuortti 1 and Nuortti 2 approx. 100 years, Kotovaara approx. 65 years and Sodankylä approx. 70-80 years. All areas had a fence that excluded all reindeer, dividing the areas to the grazed and non-grazed sites. Nuortti 1 and 2 were located close to the border of Finland and Russia and the reindeer fence along the borderline has prevented Finnish reindeer from going to Russian side for approx. 100 years. Kotovaara and Sodankylä areas were located close to the research stations and have been fenced for around 20 and 55 years, respectively. The size of the exclosures were 90 × 110 m in Kotovaara and 40 × 50 m in Sodankylä ( Supplementary Fig. S1). The distance between Nuortti 1 and Nuortti 2 areas was approx. 1 km, Nuortti areas and Kotovaara approx. 7 km, and Kotovaara and Sodankylä approx. 150 km.
The soils in the Värriö study areas were classified as haplic podzol (FAO, 1990) on sand tills, where the bulk of the mineral soil consisted of sand (Köster et al., 2017). The study sites belonged to the Pohjois-Salla reindeer herding district, with approx. 2.2 reindeer per km 2 (Turunen et al., 2016). In Sodankylä, the soil was haplic podzol on glacioflucial sandy deposit. The soil properties of each area are presented more detail in Supplementary Table S1. The main tree species in the study areas were Scots pine (Pinus sylvestris L.) and the ground vegetation consisted mainly of shrubs (mostly Vaccinium vitis-idaea L., Empetrum nigrum L. and Calluna vulgaris L.) and various mosses. Cladonia lichens were abundant on the sites where reindeer grazing was excluded (Köster et al., 2013(Köster et al., , 2015b. The vegetation biomass and tree distribution of each area are presented more detail in Supplementary  Table S2. There were no underlying permafrost in the areas. The sampling followed a split plot experiment with four different whole plots (Sodankylä, Nuortti 1, Nuortti 2 and Kotovaara) described in Köster et al. (2013Köster et al. ( , 2015bKöster et al. ( , 2017, referred here as study areas. The Sodankylä, Nuortti 1, Nuortti 2 and Kotovaara constituted of 5, 4, 6 and 4 subplots, respectively ( Supplementary Fig. S1). These subplots were divided into grazed and non-grazed sites with pre-established fences. The subplots were assigned in a systematic manner into even areas with no visible differences in stand structure, soil or topography within each whole plot. The distance between each subplot was approx. 50 m and each replicate approx. 2 m.
Soil samples were collected in June 2013 from the humus horizon, from a 0.25 m by 0.25 m quadrat area after removing the litter. For each sample, two adjacent humus samples were taken and pooled. The samples were transferred into 1.5-ml Eppendorf vials, frozen at −180°C in liquid nitrogen within a few hours after collection and transported to the laboratory on dry ice. In total, we had 38 pooled soil samples from four sampling areas, with 19 grazed and 19 non-grazed samples. The number of replicate DNA samples from Nuortti 1, Nuortti 2, Kotovaara and Sodankylä were 5, 4, 6 and 4, respectively.

Litterbag experiment and enzyme activity measurements
For the litter bag experiment, Scots pine needle litter from a pine stand from Hyytiälä (61°51′N, 24°17′E) Finland were collected in year 2012 using litter traps. The collected litter was air dried at room temperature until a constant mass was reached. The litter bags (7 cm × 10 cm) were made of 0.2 mm nylon mesh that prevented plant roots but allowed fungal hyphae to penetrate in the litter bags (Köster et al., 2015a). Approx. 5 g of dried litter were placed in the bags and the bags were buried into the organic horizon in year 2013 to the same grazed and non-grazed study areas in Värriö. Litter bags were harvested one year later and frozen after collection. Before the analysis, samples were homogenized, weighed and a subsample was taken for the enzyme activity measurements. The rest of the litter was dried in an oven at 50°C until a constant mass was reached to measure the litter dry weight and moisture content. Mass loss from the litter was determined as the difference between the initial (in 2013) and final (in 2014) dry weights of the litter.

Fungal community analysis
DNA was extracted from 0.1 g (fresh weight) of humus soil with Mobio PowerSoil DNA extraction kit (Mo Bio Laboratories, Carlsbad, CA, USA), according to the manufacturer's instructions. Bead beating was conducted with FastPrep instrument (MP Biomedicals, LLC, France) with speed 4 m s −1 for 30 s. DNA was further purified using GeneClean Turbo -kit (MP Biomedicals, LLC, France). DNA was quantified using NanoDrop 1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and the concentrations were adjusted to 10 ng µl −1 .
Fungal ITS2 region was amplified with gITS7 and ITS4 primers (Ihrmark et al., 2012), containing the appropriate 454 pyrosequencing A-and B-adapters and a unique 6-bp barcode sequences (Supplementary Table S3). The PCR was conducted in 25-µl volume, using Phusion high-fidelity DNA polymerase (Thermo Scientific, Vantaa, Finland) with 30 ng of template DNA. All samples were amplified as triplicates and the amplification was conducted as in Santalahti et al. (2016). The triplicate PCR products were pooled and purified using Agencourt AMPure XP beads (Beckman Coulter, Bernried, Germany). DNA concentrations were determined and 100 ng of each amplicon sample were combined into the sequencing library. The sequencing was conducted in the Institute of Biotechnology at the University of Helsinki using 454 GS-FLX Titanium protocol (454 Life Sciences/Roche Diagnostic, CT, USA). The sequence data are available in the European Nucleotide Archive at the European Bioinformatics Institute under study PRJEB21587 with accession numbers ERS1810502-ERS1810541.
The sequencing data were analyzed using mothur (v. 1.33.3) pipeline (Schloss et al., 2009), and the modified protocol for fungi was used as in Santalahti et al. (2016). Shortly, the sequences were trimmed and quality checked. Sequences containing any ambiguous (N) bases, homopolymers longer than eight nucleotides, or the average Phred quality score lower than 25 were discarded. All sequences were truncated to 200 bp after removing the primes and barcode sequences. Sequences were pre-clustered with 1 bp distance using a pseudo-single linkage algorithm in mothur (Huse et al., 2010) to remove sequences resulting from pyrosequencing errors. Chimeric sequences were removed using mothur-embedded UCHIME (Edgar et al., 2011). Unique sequences were pairwise aligned using the Needleman method (Needleman and Wunsch, 1970) and the obtained distance matrixes were clustered into operational taxonomic units (OTUs) using the average neighbor algorithm with 97% similarity. All global singletons were omitted due to their uncertain origin (Tedersoo et al., 2010). The sequence counts were normalized as relative abundance data within each sample by dividing the read count of each OTU by the library size obtained from the sample.
For each OTU, the representative sequences were assigned to taxa with an 80% confidence threshold using the naïve Bayesian classifier (Wang et al., 2007) and the mothur-formatted UNITE taxonomy reference database (UNITE + INSD, version 6, released 2015-03-02) in mothur. Non-fungal OTUs were removed from further analysis. For calculating the diversity estimates, the sequence data were randomly subsampled to 2500 reads per sample to ensure comparable diversity estimators across the samples (Gihring et al., 2012) using the rar-efy_even_depth function in phyloseq package (McMurdie and Holmes 2013) in R software version 3.3.3 (R Core Team 2017). The diversity estimates of observed species richness (Sobs), estimated species richness (Chao1) (Chao, 1984) and inverse Simpson index (InvSimpson, 1/ D) were further calculated using estimate_richness function.
To calculate the number of shared and unique OTUs, the OTU-data was precence/absence transformed. The percentages of the shared and unique OTUs were calculated between the grazing treatments or sampling areas, and the data were visualized using venn diagram. Fungal OTUs were further categorized into ectomycorrhizal (ECM), ericoid mycorrhizal (ERM), saprotrophic (SAP) and 'other' groups based on their functional role, by annotating the OTUs against the FUNGuild database (Nguyen et al. 2016). The group 'other' contained all other functional groups than ECM, ERM and SAP, such as endophytic, pathogenic, parasitic and lichenized fungi. Those OTUs that could not be assigned to any groups were named 'unassigned'.

Statistical analysis
Two-way ANOVA was used to analyze the statistical differences between the grazing treatments and the sampling areas for each diversity estimates (Sobs, Chao1 and Inverse Simpson index), and mass loss percentage from the needle litter. The percentage data was arcsine transformed to meet the criteria for two-way ANOVA. In the analysis, grazing was used as fixed and area as random factor. Two-way ANOVA was performed using software IBM SPSS Statistic 22 and P ≤ 0.05 was set as the limit for statistical significance.
The effect of area and grazing for the fungal communities and each fungal phylum were tested using multivariate analysis of variance (MANOVA) from the normalized sequence counts with adonis-function in vegan package (Oksanen et al., 2017) in R software (R Core Team 2017). Statistical significance of the grazing effect for the relative abundances of the most abundant genera and species, ecological groups, and extracellular enzyme activities were tested with nonparametric Kruskal-Wallis test using IBM SPSS Statistic 22.
Canonical correspondence analysis (CCA) was used to explore the relationship between the fungal community structure and driving environmental variables (all numerical variables from Supplementary Tables S1 and S2) measured in earlier studies from the same reindeer grazing sites (Köster et al. 2013(Köster et al. , 2015b using cca-function from vegan package (Oksanen et al., 2017) in R (R Core Team 2017). For the analysis, both relative abundance sequence data and environmental data were logarithmic normalized with the decostand-function and the 300 most abundant OTUs were used in the analysis. The environmental variables were chosen for the CCA by model building using ordistepfunction in vegan package. The cross-correlation between the selected variables in the final model were further tested with the vif.cca-function in vegan (variance inflation factor, VIF < 4). The statistical significance of each environmental variable and CCA-axis were tested with anova.cca function with 999 permutations. The species scores of each OTU were selected from the model using ordiselect function in goevegpackage, and those OTUs which were 80% of the most abundant and 15% of the best fit, were visualized in the CCA.
CCA was also used to explore the relationship between the extracellular enzyme activities and explanatory variables from the needle litter using vegan package (Oksanen et al., 2017). Moisture content and mass loss measured from the litter bags were used as explanatory variables. Before the analysis, the enzyme activity data were normalized by dividing each value by the total sum of the enzyme activities across all the studied samples with the decostand-function and the total-method. CCA was conducted similarly as for the sequencing data. The species scores of individual enzymes were shown and the visualized CCA was scaled according to the site scores (scaling = 1). The grouping of the variables was further visualized with ellipses representing the 75% confidence intervals using ordiellipse function from vegan package (Oksanen et al., 2017).

General information on the pyrosequencing data, species richness and diversity
The 454 pyrosequencing of the 38 samples resulted in total of 213 905 high quality fungal sequences (without singletons) that were assigned to 1 186 OTUs with 97% similarity. The grazed and nongrazed sites shared 52.8% of all OTUs, and 21.1% and 26.1% of all OTUs were unique to the sites, respectively (Supplementary Fig. S2a). The four areas Sodankylä, Nuortti 1, Nuortti 2 and Kotovaara shared 9.7% of all OTUs, and 14.8%, 7.6%, 17.5% and 6.9% OTUs were unique for the areas, respectively (Supplementary Fig. S2b). The sequencing effort of all samples are presented as rarefaction curves in Supplementary Fig. S3. Overall, reindeer grazing increased the estimated species richness (Chao1) (P ≤ 0.01, with two-way ANOVA), but had no effect on observed OTUs (Sobs) or fungal diversity (Inverse Simpson's index 1/D) (Table 1).

Litter mass loss and enzyme activity
In the litter bag experiment, the needle litter had lost 25-30% of the mass after being buried underneath the humus horizon for one year. The mass loss varied between each study area, but did not differ significantly between the grazed and non-grazed sites.
Grazing significantly affected the extracellular enzyme activities measured from the needle litter, and was the most important individual parameter explaining the variation in the data (P ≤ 0.001, F = 7.9676, anova.cca) in the CCA (Fig. 1). Grazing increased CBH I activities (P ≤ 0.05, with Kruskal-Wallis) and decreased laccase activities (P ≤ 0.05) (Supplementary Table S4). The activities of ß-xylosidase, ßglucuronidase, N-acetylglucosaminidase, ß-glucosidase and acid phosphatase did not differ significantly between the grazed and non-grazed Table 1 Diversity estimates for each study area and grazing treatment. All diversity estimates, observed (Sobs) and estimated (Chao1) species richness and diversity with the inverse Simpson index (InvSimpson) were calculated from the subsampled data (2500 sequences for each). Chao1 was statistically higher in the non-grazed compared to the grazed sites (P ≤ 0.01 with two-way ANOVA). There were no statistical differences between the grazed and non-grazed sites for Sobs and InvSimpson. SD indicates for standard deviation.

Area
Grazing n Sobs ± SD Chao1 ± SD InvSimpson ± SD  Fig. 1. CCA showing the differences in the enzyme activities between each sampling area and grazing treatment. Normalized enzyme activities were used. CCA 1 represents the constrained ordination of the data with 28.81% (P ≤ 0.001, F = 15.4715) and CCA 2 with 6.43% (P ≤ 0.05, F = 3.4535) of the total variation. The grazing effect is further visualized with ellipses representing the 75% confidence intervals. The species scores of individual enzymes are presented and the significance of each explanatory variable are calculated with anova.cca and marked with asterisk ( ** P ≤ 0.01, * P ≤ 0.05, ′P ≤ 0.1). Table S4), but the overall pattern of these enzymes were associated to the grazed sites in the CCA. Within the Nuortti 2 area, the activities of ß-xylosidase, N-acetylglucosaminidase (P ≤ 0.05, for both) and ß-glucuronidase (P ≤ 0.1) were higher in the grazed compared to the non-grazed side (Supplementary Table S4). Moisture content of the needle litter was a significant environmental parameter explaining the variation in the enzyme activity data (P ≤ 0.05, F = 4.5236, with anova.cca), (Fig. 1). Mass loss was not a significant parameter and was discarded from the analysis.

Fungal community structure
Each study area harbored unique fungal community, as the CCA showed that the fungal communities within each area clustered together (Fig. 2). Grazing was significant parameter explaining the variation in the sequence data (P ≤ 0.05, F = 1.3072, anova.cca), and the fungal communities of the grazed and non-grazed sites mainly clustered separately within each area.
In the model building for the environmental variables, ground vegetation biomass, number of small trees (lower than 1.3 m) per hectare, ground vegetation root biomass, tree root biomass, number of big trees (higher than 1.3 m) per hectare and soil temperature were significant variables explaining the data and were selected in the final model. Further analysis showed that the selected variables significantly explained the variation in the fungal community structure (Fig. 2,  Table 2). Ground vegetation biomass, number of small trees per hectare and ground vegetation root biomass were the most significant variables explaining the variation in the fungal community structure (P ≤ 0.001, anova.cca), followed by tree root biomass (P ≤ 0.01), and number of big trees, soil temperature and grazing (P ≤ 0.05), (Table 2).

Fungal community structure on phylum, genus and species level
Fungal OTUs were classified into six fungal phyla. Basidiomycota were the most abundant phylum in all grazed and non-grazed areas, with on average 60.3% of all sequences and 37% of all OTUs (Fig. 3a,  Supplementary Table S5). Ascomycota were the second most abundant phylum with 20.2% of all sequences and 30.8% of all OTUs (Fig. 3a). Ascomycota was the only phylum which had statistically significant difference for grazing (P ≤ 0.05, F = 1.9644 with MANOVA) and area (P ≤ 0.001, F = 3.7727 with MANOVA). Chytridiomycota covered 13.5% of all sequences and 2.7% of all OTUs (Fig. 3a). Former phylum of Zygomycota covered 3.7% of all sequences and 3.5% of all OTUs Fig. 2. CCA showing the differences in the fungal community structure between each sampling area and grazing treatment with environmental parameters as explanatory variables. Logarithmic normalized data from the 300 most abundant OTUs were used. CCA 1 represents the constrained ordination of the data with 9.23% (P ≤ 0.001, F = 4.0281) and CCA 3 with 5.10% (P ≤ 0.001, F = 2.2232) of the total variation of each axis. The species scores of those OTUs which were 80% of the most abundant (from the top 300 OTUs) and 15% of the best fit are presented and the significance of each environmental parameter are calculated with anova.cca and marked with asterisk ( *** P ≤ 0.001, ** P ≤ 0.01, * P ≤ 0.05). The significance of each environmental parameter are presented in Table 2.

Table 2
The significance of the environmental parameters explaining the variation in the fungal community structure (Fig. 2). The significance of each parameter are marked with asterisk ( *** P ≤ 0.001, ** P ≤ 0.01, * P ≤ 0.05).  (Supplementary Table S5). Glomeromycota and Rozellomycota were found in minority (0.02% and 0.03% of all sequences, respectively). Sequences that could not be classified into phylum level covered 2.2% of all sequences and 25.3% of all OTUs. ECM fungi dominated the humus communities in both grazed and non-grazed sites with approx. 50% of all sequences, as ERM, SAP and 'other' fungi were less abundant (4.3%, 5.3% and 0.6%, respectively) groups (Fig. 3b). Approx. 39% of all sequences could not be assigned to any functional group. The abundance of ERM fungi were higher in the non-grazed sites (P ≤ 0.05, with Kruskal-Wallis), and the abundance of other functional groups were similar in the grazed and non-grazed sites.
At genus level, Cortinarius was the most abundant genus with 24.8% of all sequences, followed by Piloderma (10.2% of all sequences), Suillus (6%), Lactarius (5%), Rhizoscyphus (3.2%), Mortierella (3%) and Hydnellum (1.8%). Grazing affected the abundance of certain fungal genera, but the differences were more obvious in the youngest exclosure areas of Kotovaara and Sodankylä (Table 3, Supplementary  Table S5). Rhizoscyphus and Oidiodendron were more abundant in the non-grazed compared to the grazed sites (P ≤ 0.05 and P ≤ 0.1, respectively with Kruskal-Wallis). Within the youngest exclusion area of Kotovaara, the abundance of Lactarius was higher and the abundance of Rhizoscyphus, Meliniomyces and Umbellopsis were lower in the grazed side compared to the non-grazed side (P ≤ 0.05, for each) (Supplementary Table S5).
Up to 24 Cortinarius species were identified, and their abundances varied between the grazed and non-grazed sites and the study areas (Supplementary Table S5). C. semisanguineus was the most abundant Cortinarius species and was present in both grazed and non-grazed sites and all study areas (Table 3, Supplementary Table S5). The second most abundant species C. caperatus was only found in the Nuortti 1 area and mainly in the non-grazed side. C. coleoptera and C. brunneus were more abundant in the grazed sites and C. pseudoduracinus was more abundant in the non-grazed sites (P ≤ 0.05, for all). C. traganus, C. testaceofolius and C. carabus were only found in the grazed sites. In the CCA, C. angustinusporus and C. neofurvolaesus associated with Sodankylä area, C. caperatus and C. gentilis associated with Nuortti 1 area, and C. suberi and C. carabus associated with Kotovaara area (Fig. 2).
The abundance of other fungal species also varied between the grazed and non-grazed sites (Supplementary Table S5). Rhizocyphus ericae (with 97% of Rhizocyphus sequences) was more abundant in the non-grazed sites (P ≤ 0.05). Lactarius rufus and L. musteus (with 82% and 15% of all Lactarius sequences, respectively) were generally more abundant in the grazed sites (P ≤ 0.1, for L. musteus). Russula decolorans, R. altorubens and R. vinosa were mainly found in the Nuortti 1 and Kotovaara non-grazed sites. In the CCA, Hydnellum aurantiacum and Tylospora asterophora associated with the Sodankylä area in the CCA (Fig. 2). L. rufus and Meliniomyces bicolor associated with the Kotovaara grazed site, and Lactarius glyciosmus associated with the non-grazed site. OTUs belonging to Mortierellales, Helotiales and Chaetohyriales associated with the Nuortti 1 area.

Discussion
To our best knowledge, our field study is the first to describe detailed fungal community structure in sites with long-term reindeer grazing and exclusion in northern boreal coniferous forests using highthroughput sequencing. Overall, reindeer grazing had significant impact on soil fungal community structure and activity of enzymes related to litter degradation, even though the decomposition rate of litter was not affected. Grazing increased the estimated species richness but had no significant effect on soil fungal diversity. Vegetation changes caused by reindeer grazing altered the fungal community structure, as multiple grazing related environmental variables significantly explained the variation of the fungal communities. On a phylum level, grazing decreased the abundance of Ascomycota and ERM fungi. Contrarily to our hypothesis, the abundance of ECM fungi was similar in the grazed and non-grazed sites and was not affected by the presence of lichens. Grazing also affected the abundance of certain fungal genera and species, but the effect was more visible in the youngest exclusion areas of Kotovaara and Sodankylä. Our study sites were located in four areas within 150 km from each other, thus our data provides a representative insight into the reindeer grazing effects in northern boreal forests in Finland.

Impact of reindeer grazing on litter decomposition and enzyme activities
Fungi produce wide range of extracellular enzymes, and are the most efficient organic matter decomposers in boreal areas. Mass loss of the needle litter was approx. 25-30%, and did not differ significantly between the grazed and non-grazed sites, after being buried in organic horizon for one year. Our results are consistent with the previous studies by Stark et al. (2010Stark et al. ( , 2002 showing that no clear grazing-related differences in litter decomposition rate were found in northern boreal forest. Contradicting result have also been obtained as reindeer grazing retarded the decomposition rate of V. myrtillus leaves in Scots pine forest in Finland (Stark et al., 2000) and fine roots in Alaskan taiga (Ruess et al., 1998), and accelerated litter decomposition rate and nutrient availability in tundra heat in northern Norway (Olofsson et al., 2004).
In our study, reindeer grazing generally enhanced the extracellular enzyme activities related to cellulose and hemicelluloses degradation and nutrient acquisition from the needle litter. Grazing significantly enhanced cellulose chain hydrolyzing CBH I activity. In addition, the activities of ß-xylosidase, ß-glucosidase, ß-glucuronidase and N-acetylglucosaminidase, as well as phosphate group cleaving acid phosphatase activities were associated to the grazed sites in the CCA. Previously, the activities of ß-glucosidase, acid phosphatase and leucine-aminopeptidase activities were higher and N-acetylglucosaminidase was lower in heavily than lightly grazed areas in a mesotrophic tundra heat (Stark and Väisänen 2014). They further suggested that grazing altered the enzyme activities through substrate availability for soil microorganisms (Stark and Väisänen 2014). In our Table 3 The relative abundance of the most abundant genera and species in the grazed and non-grazed sites. Statistical significance of the grazing effect are calculated with non-parametric Kruskal-Wallis test and marked with asterisk ( * P ≤ 0.05, • P ≤ 0.1). SD indicates for standard deviation, n = 19.

Genus
Grazed ± SD Non-grazed ± SD study, moisture content of the needle litter was also significant variable explaining the variation in the enzyme activities, which was in line with previous studies where soil moisture content was found to correlate positively with enzyme activities (Baldrian et al. 2010a,b). Contrarily to other measured enzyme activities, grazing significantly decreased the activity of lignin modifying laccase. Previously, grazing was found to decrease soil microbial respiration and activity, which was suggested to result in decreased decomposition rate (Väre et al. 1996, Stark et al. 2003. Increased N availability from reindeer feces in heavily grazed areas have also been found to reduce the decomposition rate of recalcitrant organic matter (Craine et al 2007), microbial biomass and activity of ligninolytic enzymes (Sinsabaugh 2010;Männistö et al. 2016). Although we did not measure N availability, microbial biomass was generally lower in the grazed areas (Köster et al., 2015b). Further, as several lichens (including some Cladonia species) have shown to possess laccase activities (Beckett et al. 2014;Zavarzina and Zavarzin 2006), we cannot rule out the possibility that the higher abundance of lichens in the non-grazed sites could have affected to the higher laccase activities on these sites, although the laccase activities from non-Peltigeralean lichens have generally been low (Beckett et al. 2013).
There are numerous contradictory results concerning the effects or reindeer grazing on litter and SOM decomposition rate, activity of various extracellular enzymes and soil microbial activity. It is evident that several factors are simultaneously affecting soil processes and therefore, the grazing effect should be investigated more thoroughly in a multidisciplinary manner. It is also noted that our one-year litter bag experiment might not had enough time to induce clear differences in the litter decomposition rate, even though there were clear differences in the enzyme activities between the grazed and non-grazed sites. As grazing decreases the number of young tree seedlings, it may affect the litter input, and consequently alter litter and SOM quality. Nevertheless, these processes take time as the tree growth and SOM decomposition processes in northern boreal forests are extremely slow.

Impact of reindeer grazing on soil fungal community structure
Reindeer grazing altered the fungal community structure, as grazing was significant parameter explaining the variation in the CCA. However, area was more important parameter, as the communities clustered together based on the study areas in the CCA, and all study areas harbored a unique species composition. These findings are in line with earlier studies showing that spatial heterogeneity is an important factor affecting microbial communities more than e.g. grazing (Stark et al., 2008;Sørensen et al., 2009).
The relationship between fungal community structure and driving environmental variables showed that many of the studied variables were significant and explained the variation in the fungal community structure. Previously, reindeer grazing was found to significantly reduce lichen and total ground vegetation biomasses (mainly due to the decrease in lichen biomass) (Köster et al., 2013(Köster et al., , 2015b. Contrarily to our hypothesis, the abundance of ECM fungi in general varied inconsistently between the grazed and non-grazed sites, and the presence of lichens did not suppress most ECM fungi. The vegetation changes due to reindeer grazing were not only limited to the reduction of lichen biomass, but also to the reduction of small tree seedlings (Köster et al., 2013(Köster et al., , 2015b. The inconsistency in the number of big trees between the treatments may have caused differences in the abundance of root and fine root biomass, and thus may also reflect to the abundance of ECM fungi. Even though the fine roots were not quantified in this study, the tree root biomass significantly explained the variation in the fungal community structure. Therefore, further studies are needed to validate this via quantifying the fine root biomass.

Genus and species level changes
Cortinarius was the most abundant genus by the number of sequences and OTUs, and the abundance of Cortinarius species varied between the grazed and non-grazed sites. Genus Cortinarius are abundant and species rich ECM fungi in boreal forest soils (Bödeker et al., 2014;Lindahl et al., 2010Lindahl et al., , 2007Santalahti et al., 2016;Sun et al., 2015), in Arctic tundra (Deslippe et al., 2011) and in eutrophic heaths and meadows in Fennoscandia (Lye, 1975). Cortinarius are strongly proteolytic fungi that are adapted to the typical N-limited conditions of boreal forests and derive C from their host for biomass and extracellular enzyme production (Deslippe et al., 2011). In the Scandies Mountain Range, Cortinarius was found to be more abundant in non-grazed exclosures than grazed control plots (Vowles et al., 2018).
According to our previous study from Värriö, genus Cortinarius was clearly associated with the functional gene profile observed at the old growth forest site , suggesting that Cortinarius also are highly important genus for the functioning of the old growth forest soils in the area, close to our study sites. Furthermore, ECM Cortinarius species have been detected to possess oxidative enzymes, such as laccase and manganese peroxidase (MnP) (Bödeker et al., 2014;Courty et al., 2006), and it has been suggested to have significant role in the decomposition of complex SOM. Several Cortinarius species have been demonstrated to co-localize with high MnP activity (Bödeker et al., 2014). In our study, some species e.g. C. caperatus, C. biformis and C. suberi were in general more abundant in the non-grazed sites with higher laccase activity. In addition to genus Cortinarius, several other ECM fungi have also been demonstrated to produce various oxidative enzymes, including genera Piloderma, Lactarius and Suillus (Bödeker et al., 2014;Heinonsalo et al., 2012;Kohler et al., 2015;Shah et al., 2015), which all were abundant in our study sites. Since fungi possess large pool of various lignin modifying enzymes, further study to detect grazing induced differences in the patterns to modify lignin to its subunits would be needed.
Other highly abundant genera in our study were Piloderma, Suillus, Lactarius, and Rhizoscyphus, which all are typical mycorrhizal fungi of boreal forest soils (Markkola et al., 2002;Santalahti et al., 2016;Sun et al., 2015). Genus Rhizoscyphus and species R. ericae, which forms ERM with roots of Ericaceae plants (Smith and Read, 2008), were statistically more abundant in the non-grazed sites, even though grazing did not affect shrub biomass (Köster et al., 2013(Köster et al., , 2015b. It is also possible that the roots and especially the fine roots of Ericaceae and other shrubs are more susceptible for the trampling of reindeer than e.g. tree roots that mainly are located deeper in the soil (Makkonen and Helmisaari, 1998;Helmisaari et al., 2007).
Suillus variegatus, which was the most abundant species among genus Suillus in our study areas, and Piloderma croceum has previously been reported to suffer from lichen extracts (Brown and Mikola, 1974;Markkola et al., 2002), but also stimulated by reindeer lichen Cladonia rangiferina (Brown and Mikola, 1974). In our study, genus Suillus and Piloderma varied inconsistently between the grazed and non-grazed sites and did not seem to suffer from lichens. In contrast, the ECM genus Lactarius and the most abundant species L. rufus and L. musteus were on average more abundant on the grazed sites and might have suffered from the lichens.

Conclusions
Our study is the first to describe detailed fungal community structure using a 454-pyrosequencing analysis in reindeer grazed and nongrazed sites in northern boreal coniferous forest soils. Reindeer grazing altered significantly the soil fungal community structures and the abundance of Ascomycota, as well as certain fungal genera and species. Furthermore, the activities of extracellular enzymes related to litter degradation were also affected by grazing, even though the measured litter mass loss did not differ between the grazing treatments after one year. Our data indicates that with longer time scales, grazing may affect litter decomposition through changes in fungal community structure and enzyme activities in the northern boreal forest soils.