Reindeer control over subarctic treeline alters soil fungal communities with potential consequences for soil carbon storage

The climate‐driven encroachment of shrubs into the Arctic is accompanied by shifts in soil fungal communities that could contribute to a net release of carbon from tundra soils. At the same time, arctic grazers are known to prevent the establishment of deciduous shrubs and, under certain conditions, promote the dominance of evergreen shrubs. As these different vegetation types associate with contrasting fungal communities, the belowground consequences of climate change could vary among grazing regimes. Yet, at present, the impact of grazing on soil fungal communities and their links to soil carbon have remained speculative. Here we tested how soil fungal community composition, diversity and function depend on tree vicinity and long‐term reindeer grazing regime and assessed how the fungal communities relate to organic soil carbon stocks in an alpine treeline ecotone in Northern Scandinavia. We determined soil carbon stocks and characterized soil fungal communities directly underneath and >3 m away from mountain birches (Betula pubescens ssp. czerepanovii) in two adjacent 55‐year‐old grazing regimes with or without summer grazing by reindeer (Rangifer tarandus). We show that the area exposed to year‐round grazing dominated by evergreen dwarf shrubs had higher soil C:N ratio, higher fungal abundance and lower fungal diversity compared with the area with only winter grazing and higher abundance of mountain birch. Although soil carbon stocks did not differ between the grazing regimes, stocks were positively associated with root‐associated ascomycetes, typical to the year‐round grazing regime, and negatively associated with free‐living saprotrophs, typical to the winter grazing regime. These findings suggest that when grazers promote dominance of evergreen dwarf shrubs, they induce shifts in soil fungal communities that increase soil carbon sequestration in the long term. Thus, to predict climate‐driven changes in soil carbon, grazer‐induced shifts in vegetation and soil fungal communities need to be accounted for.


| INTRODUC TI ON
Arctic and subarctic regions are rapidly warming, which is leading to northward and altitudinal range expansion of shrubs and trees (Myers-Smith et al., 2019;Tape et al., 2012;Terskaia et al., 2020).
Although this 'shrubification' likely contributes to higher primary production (Zhang et al., 2014), it could also accelerate fungal-driven soil organic matter (OM) decomposition and result in a large net loss of ecosystem carbon and exacerbate climate change (Hartley et al., 2012;Parker et al., 2015). At the same time, virtually all of the Arctic experiences various degrees of vertebrate and invertebrate herbivory (Barrio et al., 2016;Rheubottom et al., 2019), which is known to hamper tree and shrub establishment and growth (Dahlgren et al., 2009;Kumpula et al., 2011;Sundqvist et al., 2019), induce tree mortality (Jepsen et al., 2013) and keep tree-lines lower than their altitudinal limits (Cairns & Moen, 2004;Oksanen et al., 1995;Van Bogaert et al., 2011). Thus, to truly predict the realized vegetation ranges and ecosystem processes in warmer temperatures, the past, current and future herbivore densities need to be accounted for (Schmitz et al., 2018).
Reindeer (Rangifer tarandus L.), also known as the caribou in North America, is the most wide-spread vertebrate herbivore in the Arctic (Bernes et al., 2015) with wild or semi-domesticated reindeer herds roaming over the vast majority of Arctic tundra (Conservation of Arctic Flora and Fauna (CAFF), 2001). These animals adapt to the pronounced seasonality of the Arctic by altering their metabolic efficiency between seasons and by migrating large distances to optimize forage intake according to availability (Bernes et al., 2015). During winters, reindeer lichens (Cladonia) are a central nutrition source for most reindeer herds, whereas during summers, reindeer browse on the leaves of shrubs and deciduous trees and graze on graminoids (grasses, sedges and rushes) and herbs.
Here, we refer to the combined effects of both foraging behaviours as grazing. Depending on the season and reindeer abundance, reindeer may alter tundra vegetation and, in certain conditions, even shift the ecosystem to an alternative state with altered process rates and limited capacity to recover to the earlier state (Egelkraut et al., 2018;van der Wal, 2006;Ylänne et al., 2020). Observations of reindeer-induced state shifts have led to suggestions that reindeer management (Bråthen et al., 2017) and even reintroductions of large grazers (Olofsson & Post, 2018) could be used to reduce some of the unwanted effects of a warmer climate, particularly the shrubification. However, the consequences of such interventions on soil carbon storage remain largely uncertain (Olofsson & Post, 2018) and are likely dependent on the characteristics of the ecosystem, such as the dominant plant traits (sensu Wardle et al., 2004).
For example, ecosystems dominated by fast-growing plants with high nitrogen content are likely to benefit from the grazing-induced high return of labile fecal material to the soil, leading to herbivory promoting faster soil nutrient turnover and higher productivity rates (Wardle et al., 2004). In contrast, ecosystems dominated by slow-growing plants could respond to herbivory by increasing the amount of defence compounds, which would reduce litter quality and ultimately lead to reduced soil nutrient availability and productivity in grazed systems (Wardle et al., 2004). Field surveys from the tundra show evidence of reindeer-induced shifts in both directions, with either positive (Egelkraut et al., 2018;Olofsson et al., 2004;Stark & Väisänen, 2014) or negative (Sundqvist et al., 2019;Vowles et al., 2017;Ylänne et al., 2015) consequences on plant palatability and nitrogen content. Which one of these alternative scenarios takes place is likely to determine grazing impacts on soil microbial communities and the ecosystem processes they regulate, such as soil carbon sequestration (Wardle et al., 2004).
Grazers' control over soil carbon sequestration could be particularly important in situations where grazing favours the dominance of evergreen shrubs over deciduous trees and shrubs. This could specifically be relevant in Fennoscandia where the deciduous genus Betula has been suggested to be key for ecosystem carbon storage (Hartley et al., 2012;Parker et al., 2018). Although the high-stature Betula species store large amounts of carbon in their biomass, total ecosystem carbon storage is lower in both mountain birch (Betula pubescens ssp. czerepanovii) forests and dwarf birch (B. nana) shrublands compared with evergreen shrub-dominated tundra heaths due to lower soil carbon stocks in Betula-dominated ecosystems (Hartley et al., 2012;Parker et al., 2018). The differences in soil carbon have been suggested to be linked to the mycorrhizal fungi associated with the dominant vegetation (Parker et al., 2015;Clemmensen et al., 2021). Although ectomycorrhizal fungi (EcM) associated with birch may generally grow and turn over biomass fast, many of the ericoid mycorrhizal fungi (ErM) associated with evergreen shrubs have melanised cell walls that are relatively resistant to decomposition (Clemmensen et al., 2015;Fernandez & Kennedy, 2018;Read & Perez-Moreno, 2003). In addition, certain EcM fungi are known to efficiently access and mine nitrogen from OM by using oxidative enzymes (Bödeker et al., 2014;Kyaschenko et al., 2017) and by traits adapted to long-distance nutrient transport, such as cords (Clemmensen et al., 2015). This contributes to EcMdominated ecosystems generally maintaining higher rates of OM decomposition-and lower rates of humus accumulation-compared with ecosystems dominated by ErM fungi (Clemmensen et al., 2015. The difference between EcM-dominated mountain birch forests and ErM-dominated tundra heaths could be particularly strong if the difference in dominant mycorrhizal association is accompanied by a shift in free-living saprotrophic activity . Due to the higher litter quality of deciduous versus evergreen leaves (Cornwell et al., 2008) and the soil community adapted to this K E Y W O R D S Arctic shrubification, Betula pubescens ssp. czerepanovii, fungal community, grazing, ITS2, Rangifer tarandus, subarctic tundra, tree-line (Parker et al., 2018), the activity of saprotrophic fungi is generally higher in mountain birch forests compared with tundra heaths.
Given what is known about the role of herbivores in climatedriven vegetation changes and their potential to promote an evergreen shrub-dominated plant community, herbivores' potential to increase tundra soil carbon stocks via promoting the dominance of ErM has been proposed ), yet evidence for this mechanism is lacking. Here, we test this idea across two adjacent grazing regimes, which have experienced either only winter grazing or year-round grazing by reindeer due to long-standing differences in Finnish and Norwegian reindeer husbandry. The winter-only grazing regime is characterized by a low-stature mountain birch woodland, whereas summer browsing in the year-round grazing regime (YGR) has hindered tree establishment (Biuw et al., 2014;Kumpula et al., 2011) turning this area into a savannah-like open tundra heath with individual mountain birches scattered very sparsely in the landscape (Oksanen et al., 1995). Here, we firstly investigate the impacts of grazing regime and mountain birch vicinity on the abundance, diversity and community composition of the soil fungal community and secondly explore how the soil fungal communities relate to vegetation and, ultimately to soil carbon stocks. Specifically, we test the hypotheses, that (H1) year-round grazing by reindeer and greater distance to mountain birch trees increase the relative abundance of ErM, while decreasing EcM and saprotrophic fungi, and that (H2) ErM-dominated communities are associated with higher organic soil carbon stocks than EcM-dominated communities, particularly if (H3) the EcM fungi belong to the cord-forming exploration type and/or if (H4) the mycorrhizal shift from EcM to ErM is accompanied with lower abundance of free-living saprotrophs.

| Study site, experimental design and soil sampling
The study was carried out at the tree-line ecotone in the area of www.senor ge.no). The landscape consists of mountain birch forests, shrublands and tundra heaths on uplands and fens in depressions (Kitti et al., 2009). In the late 1950's, a border fence was built to separate the reindeer herds of the two countries (Forbes & Kumpula, 2009). Due to the migratory herding tradition in Norway, where reindeer move between coastal summer grazing grounds and inland winter grazing grounds, the Norwegian side of the Jávrrešduoddarat area is grazed only during snow-covered winter months. In contrast, the Finnish side of the Jávrrešduoddarat area is grazed all year round with the grazing pressure in summer being approximately 10-12 reindeer per square kilometre (Kitti et al., 2009). As reindeer exert less control over vegetation during winter months, when soil and vegetation are mostly protected by snow and ice (Kumpula et al., 2011), the year-round grazed side is more exposed to reindeer browsing, trampling and fertilization by excreta than the winter grazed side (Forbes & Kumpula, 2009). This long-term grazing difference is reflected in higher lichen cover in the winter grazing regime (WGR) compared with the YGR (Cohen et al., 2013;Forbes & Kumpula, 2009) and in the abundance of mountain birch (B. pubescens ssp. czerepanovii (Orlova) Hämet-Ahti) at the treeline ecotone (Oksanen et al., 1995). Although ~70 years of summer browsing by reindeer in the YGR has hampered regeneration of mountain birch stands (Biuw et al., 2014) and resulted in scattered trees occurring sparsely in a savannah-like open landscape, the WGR is characterized by high abundance of mountain birch (Oksanen et al., 1995).
Here, we compared the open woodland in the YGR to the denser mountain birch woodlands in the WGR by sampling three 350 × 200 m sized blocks, separated by ~1 km, along the border fence. All blocks had one side in the YGR and one side in the WGR, with no difference in topography, altitude or slope between the two sides of the fence. The average densities of mountain birch were 213 ± 115 and 3 ± 2 trees ha −1 in the WGR and YGR, respectively, with slightly taller (3.3 ± 0.5 m), one-or two-stemmed trees in the YGR and shorter (2.6 ± 0.8 m), multi-stemmed trees in the WGR (see Figure S1). In both grazing regimes, the field layer was dominated by evergreen and deciduous dwarf shrubs, such as Empetrum nigrum L.  Figure   S2a). At each coordinate, we collected soil samples below the canopy of the nearest mountain birch and the nearest point that was at least three meters away from the nearest mountain birch ( Figure   S2b), where the low-stature shrubs dominated the vegetation. This study design amounted to a total of 72 samples (3 blocks × 2 grazing regimes × 6 coordinates × 2 distances to mountain birch trees).
As this study is based on a singular border fence that separates the grazing regimes of the neighbouring countries, a risk of pseudoreplication affecting the results exists. Yet, we have minimized potential dependencies caused by this design by keeping a higher distance between consecutive plots on each side of the fence than between the plot pairs across the fence (Olofsson et al., 2004).
Soil samples were collected on the 28th and 29th of September 2019 by taking three to four soil cores (Ø 5 cm) from the upper 5 cm of the soil. The sample included the litter and humus layers in each plot, and both live vegetation and mineral soil was carefully removed, and the depth of each sample was measured. While sampling, hummocks were avoided. As hummocks typically form around tree-stems, the soil cores close to mountain birches were collected below the outer border of the canopy. Within 48 h of collection, soil samples were weighed and divided into two equally sized halves and stored at −18℃. One half of the soil sample was used to determine soil water (SWC, drying at 105℃ for >18 h) and OM content (loss on ignition at 405℃, 5 h). The other half was freeze-dried (LYO GT2; LYOVACTM; GEA Group) and mixed manually during which large roots and stones were removed. A subset (2 ml) of the mixed soil was homogenized with a ball mill (MM200; Retsch) before DNA extraction and chemical analyses. Soil carbon and nitrogen contents were analysed with vario Max CN (Elementar) and used to calculate organic soil C:N ratio and C and N stocks (SOC and SON, kg m −2 ).

| DNA extraction and quantitative PCR (qPCR)
Total DNA was extracted from 100 mg of dried soil, using the NucleoSpin® Soil kit (Macherey-Nagel) and the DNA concentration was estimated spectrophotometrically (NanoDrop 2000c; Thermo Scientific) and diluted 1:20. Polymerase chain reaction (PCR) inhibition was tested by spiking each sample with a known amount of the pGEM-16S plasmid (Promega) and amplifying a region on the plasmid using the primers M13F and M13R on a CFX Connect Real-Time System (Bio-Rad). No inhibition was detected in the samples. Copy numbers of the fungal ITS2 region were then quantified in duplicates of all samples using the fungal primers gITS7 (Ihrmark et al., 2012), ITS4 (White et al., 1990) and ITS4arch (Sterkenburg et al., 2018) on the same Real-Time System. The reaction mix (20 μl) consisted of 0.5 μM gITS7, 0.3 μM ITS4, 0.15 μM ITS4arch, 0.1% bovine serum albumin, 1 × iQ-SYBR Green (Bio-Rad) and 2 μl of diluted DNA template. The PCR cycling conditions were an initial denaturation at 95℃ for 5 min, followed by 40 cycles of denaturation at 95℃ for 15 s, annealing at 56℃ for 30 s, synthesis at 72℃ for 40 s and at 78℃ for 5 s. Standard curves were based on serial dilutions of linearized plasmids with cloned ITS2 fragments from Pilidium concavum. Total ITS2 copy numbers were corrected for the proportion of non-fungal (mainly plant) reads based on sequencing of the same region (see below).

| Soil fungal community
The fungal ITS2 region was PCR-amplified in duplicates from each sample using the same primers as above. Length distribution of the mix was assessed using a BioAnalyzer 2100 (Agilent Technologies) with a 7500 DNA chip. Samples were sequenced at SciLifeLab NGI with PacBio Sequel I (Pacific Biosciences) using one single-molecule real-time cell, after addition of sequencing adapters by ligation.

| Bioinformatic analyses
The raw sequence data, containing 372,213 reads, were qualityfiltered and clustered using the SCATA pipeline (https://scata.mykop at.slu.se). ITS2 sequences of <200 bases in length, with an average quality score of less than 20 or containing individual bases with a quality score of less than 7 were removed. Sequences were screened for primers and identification tags and reverse complemented if necessary, and sequences with less than 90% primer or 100% identification tag match were removed. In total, 223,814 sequences passed the quality control. After the collapse of homopolymers to 3 bp, sequences were clustered into species-level operational taxonomical units (OTUs) using single-linkage clustering (U SEARCH ; Edgar, 2010) with a threshold distance of 1.5% , mismatch penalty of 1, gap open penalty of 0, gap extension penalty of 1 and end gap penalty of 0 in the pair-wise comparisons. Genotypes occurring only once in the whole data set were also removed along with sequences with non-matching sample tags. In total, we obtained 160,683 high-quality, non-unique reads which clustered into 992 OTUs.
For taxonomic identification, reference sequences from the UNITE database (Abarenkov et al., 2010) and relevant reference data sets publicly available in SCATA were included in the clustering in SCATA. These provided identifications based on the same criteria as used for the clustering, without affecting it. In addition, the most common genotype of each OTU was subjected to taxonomic assignment in a PROTAX analysis with 50% probability in the PlutoF platform (Abarenkov et al., 2018). Based on taxonomic affiliation and match to references associated with particular substrates (Clemmensen et al., 2013), OTUs were separated into the following functional guilds: rootassociated ascomycetes, root-associated basidiomycetes, moulds, yeasts, other saprotrophs (including litter-associated fungi) and lichenized fungi. The root-associated fungi were further classified into EcM (including all fungal OTUs forming ectomycorrhizal symbiosis) and ErM (excluding the OTUs that also form EcM). EcM fungi were further subdivided according to exploration type (Agerer, 2006) to cordformers and simple mycelia as defined by Clemmensen et al. (2015).
After removal of non-fungal OTUs, 106,872 reads remained. To reduce the risk of random errors induced by multi-primer artefacts (Tedersoo et al., 2018), all observations with four reads or less were omitted; subsequently, four reads were subtracted from each of the remaining OTUs to restore proportionality of the rarest species in relation to the rest of the community. This resulted in a final data set with a total of 82,845 reads and 364 OTUs, and an average of 1167 (222-8106) reads per sample, used in subsequent statistical analyses. In total, 61% of the OTUs could be assigned to a functional guild, representing 86% of all reads.

| Estimated plant abundance
A vegetation survey was conducted on the study area, but it did not evaluate tree abundance and ground layer cover in the same plots as used for soil sample collection. Therefore, we used plant reads captured by the ITS2 sequencing as a proxy of plant root abundance in the sampled soils (see a comparison of vegetation survey and DNA-based data in the Method S1 description). The plant reads were converted to reads per g DW soil by multiplying their relative abundance by the total ITS2 copy numbers (qPCR) in each sample, and classified into the following groups: Betula spp., Vaccinium spp., other ericaceous dwarf shrubs (E. hermaphroditum, P. caerulea and C. vulgaris) and bryophytes (Polytrichum and Dicranum spp.). These four groups contributed 89.0 ± 6.1% of plant cover in the blocks (see Method S1).

| Statistical analysis
All statistical analyses were conducted using R version 3.6.3 (R Core Team, 2020) using the package ggplot2 ( To test the impacts of grazing regime, tree vicinity and their interaction on the relative abundance of taxonomic groups (phyla, subphyla, classes, orders and species), functional guilds and orders and species within functional guilds, multivariate generalized linear models (package mvabund v4.1.3; Wang et al., 2012) were built with block set as strata. Non-rarefied data was used for calculating relative abundances, these were rounded up to closest 0.01% and only subgroups with a relative abundance of more than 0.5% were included in the tests. Models were fitted with a negative binomial distribution using 999 bootstrap iterations and reported with log-likelihood ratio statistics. To test the impacts of grazing regime, tree vicinity and their interaction on finer functional groupings, the previous models were repeated by first replacing the root-associated guilds with the two mycorrhizal types, and subsequently by replacing EcM fungi with the two exploration types. In all functional guild models, the model was first run including the unknown function or type, and subsequently without the unknown groups (Table S1). As the initial test on taxonomy and functions showed strongest responses on order and species levels with no added impact of mycorrhizal and exploration type (Table S1), only orders and species were tested within the five functional guilds. To reveal grazing and tree impacts on the sub-groups, the multivariate generalized linear models were followed by univariate tests with unadjusted p-values (summarized for phyla, order and species in Table S2, for functional guilds in Table S3, and for orders and species within the functional guilds in Table S4).
Impacts of grazing regime and tree vicinity on fungal community composition were further visualised using non-metric multidimen-  Table S5).
To better understand the connections between the plant and fungal communities, we assessed how the captured abundance of Betula sp., Vaccinium sp., other ericaceous species and bryophytes correlated with fungal guilds, mycorrhizal types and EcM exploration types (Table   S6). For the analysis, all the saprotrophic guilds, that is, moulds, yeasts and other saprotrophs, were merged into one common group of saprotrophic fungi. All correlations included Hellinger-transformed relative abundances of the fungal groups (out of total fungal reads) and square root-transformed abundances of plant reads per DW soil.
To test whether fungal abundance and community composition explained soil organic carbon (SOC) stocks, linear models of SOC with fungal abundance (log-transformed ITS2 copy numbers) and fungal functional groups (Hellinger-transformed) as explanatory factors were built. The saprotrophic guilds, that is, moulds, yeasts and other saprotrophs, were merged for this analysis. To shed light on the taxonomical changes behind the functional guild shifts, we further inspected correlation among SOC and fungal order and species, limited to groups with relative abundance > 0.1. With all response variables, we first tested how the abundances, and the fungal functional and taxonomic groups explained SOC (Table S7) Table S8.

| Vegetation and soil characteristics in the study area
Based on the captured plant reads in the soil samples, the YGR was characterised by higher abundance of evergreen ericaceous dwarf shrubs, that is, E. hermaphroditum, P. caerulea and C. vulgaris, and bryophytes compared with the WGR (Table 1). The soil organic layer in both grazing regimes was on average 3 cm thick, and soil bulk density, water content and OM content did not differ between the grazing regimes (Table 1; Figure   S3). Yet, soil C:N ratio was higher in the YGR (average ± SD 27.8 ± 2.8) compared with the WGR (25.3 ± 3.3; F 1,65 = 11.63, p = 0.001). Regardless of grazing regime, soil samples collected in the vicinity of mountain birch were characterised by higher abundance of Betula and Vaccinium reads and higher soil OM content ( Figure S3; Table 1).

| Fungal abundance, diversity and community composition
Fungal abundance in the organic soil layer was higher in the area exposed to year-round grazing (Figure 1a; Table 1). However, the YGR had lower species richness and evenness ( Figure S4b,c; Table 1), resulting in a decreased effective diversity (Figure 1b; Table 1). All diversity indices were indicatively lowest in the YGR >3 m away from birches (Figure 1b; Figure S4b-d), yet the interaction between grazing regime and tree vicinity was significant only for the inverse Simpson diversity index that was higher away from the birches in the WGR than away from the birches in the YGR (Table 1; Figure S4d).
The grazing regime also affected the taxonomic and functional composition of the fungal community, with grazing impacts evident on phylum, order and species levels ( Table S1). The orders Mortierellales, Mucorales, Leucosporidiales, Thelephorales and Eurotiales were significantly or close-to-significantly more abundant in the WGR (Table S2; Figure 2). Instead, the phylum Asomycota, the orders Archaeorhizomycetales and Lecanorales and the most abundant species Pezoloma ericae were more abundant in the YGR (Figure 2; Figure S5; Table S2). Tree vicinity had little effect on the taxonomic composition, yet tree vicinity interacted with grazing to explain relative abundances of fungal orders and species (Tables S1 and S2; Figure 2).
For fungal functional guilds, the grazing difference was most evident in moulds, yeasts and lichenized fungi (Figure 3; Table S3). The relative abundances of moulds and yeasts were significantly higher The samples collected >3 m from mountain birches differ between grazing regimes.

TA B L E 1 Linear mixed-model results
on the impact of summer grazing, tree vicinity and their interaction on the organic soil layer, soil C and N stocks, vegetation abundance, fungal abundance and diversity. Significant effects are bolded and arrows indicate increase (↑) or decrease (↓) towards the year-round grazing regime or closer vicinity to mountain birch. In case of a significant interaction between grazing regime and tree vicinity; a least square means post hoc test was performed, and the pair-wise differences are found in the footnotes in the WGR than in the YGR (F 1,69 = 8.50, p = 0.009 for moulds and F 1,69 = 5.76, p = 0.017 for yeasts), whereas lichenized fungi were more abundant in the YGR (F 1,69 = 13.40, p = 0.002; Figure 3).
Furthermore, the abundance of lichenized fungi increased away from the trees in the WGR but decreased away from trees in the YGR (Figure 3; Table S3). The abundances of root-associated guilds were not significantly different between the grazing regimes nor dependant on the tree vicinity ( Figure 3; Table S3). Abundances of ErM and EcM fungi strongly reflected the abundances of root-associated ascomycetes and basidiomycetes, respectively ( Figure S6; Figure 3).
Most EcM fungi at our site had simple mycelia, with cord-forming EcM contributing less than 1% of all the fungal reads ( Figure S6). The inclusion of mycorrhizal or exploration types did not significantly increase the explanatory power of grazing or tree vicinity in the generalized linear model of the functional guild composition (Table S1).
In the NMDS, grazing exerted the strongest impact on the community composition (F 1,67 = 2.99, p = 0.001; Figure 4). Fungal communities in the YGR aligned with positive NMDS2 axis values in parallel to higher soil C:N ratios, lower OTU richness and diversity (Figure 4a,b) and higher abundance of ErM fungi (Figure 4c). The fungal community in the WGR was characterized by moulds and yeasts (Figure 4a,b), particularly the mould species within Mortierella (from the order Mortierellales), Umbelopsis (Mucorales) and Penicillium (Eurotiales) and the yeast orders Leucoporidiales and Saccharomycetales ( Figure S8b; Tables S4 and S5). No impact of grazing or tree vicinity was found on other saprotrophic fungi ( Figure 3; Table S3).
Tree vicinity also affected fungal community composition in the NMDS (F 1,67 = 1.62, p = 0.020), although this effect was partly masked by its interaction with the grazing regime (F 1,67 = 1.29, p = 0.091) indicating that tree vicinity was more important under the YGR (Figure 4a). Fungal communities away from birches in the year-round grazed area were typically dominated by rootassociated ascomycetes (Figure 4b), particularly belonging to the most abundant orders, Helotiales and Archaeorhizomycetales ( Figure S8a; Table S5). Additionally, fungal abundance F I G U R E 2 Relative abundance of the most common fungal orders (relative abundance > 0.5%; a), and the response of each order to mountain birch vicinity (filled bars = in close vicinity, and open bars >3 m from mountain birch) and grazing regime (b). Orders are colourcoded according to Phylum with yellows = Mucoromycota, greens = Basidiomycota and reds = Ascomycota and with light grey parts in (a) indicating the pooled relative abundances of the less common Basidiomycota and Ascomycota orders, in respective order. Bars represent mean ± 95% confidence interval in (b). Significant and close-to-significant differences among grazing regimes and tree vicinities are highlighted above the panels (ʹp < 0.1, *p < 0.05, **p < 0.01, ***p < 0.001; see Table S2 for further statistical parameters). Archaeorhiz., Archaeorhizomycetales; Chaetothy, Chaetohyriales; Lecan. und., undefined order in the class Lecanoromycetes; Leotiomyc. und., undefined order in the class Leotiomycetes; Leucosp., Leucosporidiales; WGR, winter grazing regime; YGR, year-round grazing regime [Colour figure can be viewed at wileyonlinelibrary.com] (ITS2 copies m −2 ) increased away from trees ( Figure 4b; Table S5).

F I G U R E 1 Fungal abundance (a) and effective diversity (b) in the organic soil layer in the vicinity (black bars) and >3 m away from mountain birch trees (white bars) in the winter grazing regime (WGR) and year-round (YGR) grazing regime. Bars present mean and 95% confidence interval
Fungal communities in the vicinity of trees instead had higher abundances of root-associated basidiomycetes (Figure 4b), particularly the EcM orders Russulales and Cantharellales (Figure 2; Figure S8a; Table S5).

| Associations between fungal community and vegetation DNA
Across all samples, the abundance of Betula reads in the soil correlated positively with the relative abundance of guilds associated with deciduous shrubs and trees (i.e. root-associated basidiomycetes, EcM and the two EcM exploration types, Figure S7b,c; Table   S6). Particularly, more EcM were found in plots with more Betula reads (F 1,69 = 6.67, p = 0.012, r 2 = 0.09; Figure S7c; Table S6). Both ErM-associated plant groups, Vaccinium sp. and other ericaceous shrubs, correlated positively with root-associated ascomycetes, while correlating negatively with root-associated basidiomycetes, EcM and saprotrophic fungi ( Figure S7e,h; Table S6). The abundance of lichenized fungi was positively linked to bryophyte abundance (F 1,69 = 10.65, p = 0.002, r 2 = 0.13), whereas saprotrophic fungi correlated negatively with bryophyte abundance ( Figure   S7k; Table S6).

F I G U R E 3
Relative abundances of fungal guilds (mean ± 95% confidence interval) in close vicinity (filled bars) and >3 m away from mountain birches (open bars) in the winter grazing regime (WGR) and year-round (YGR) grazing regimes (a), and the relative abundance of these guilds across all treatments (b). Significant differences among grazing regimes and tree vicinities are highlighted above the panels (* p < 0.05, ** p < 0.01, *** p < 0.001; see Table S3 for further statistical parameters). Other saprotr., other saprotrophs; Root-ass. ascom., root-associated ascomycetes; Root-ass. basidiom.

| SOC and its relation to fungal abundance and community composition
Organic soil stored on average (±SD) 3.18 ± 1.73 kg C m −2 in the WGR and 3.92 ± 2.18 kg C m −2 in the YGR (Figure 5a), yet the difference between the grazing regimes was not significant (F 1,65 = 2.55, p = 0.115). Instead, closer vicinity to trees tended to increase soil carbon stocks (F 1,65 = 3.15, p = 0.081) aligning with the observed pattern in soil nitrogen stocks ( Figure S9; Table 1). In the NMDS, SOC increased in the same direction as fungal abundance (Figure 4a,b), and these variables correlated also in the linear model (Figure 5b; Table   S7). Furthermore, the same taxonomic groups that aligned with fungal abundance along the NMDS2-axis ( Figure S8) related to soil carbon storage (Figure 5c; Table S7). High abundance of saprotrophic fungi, that is, moulds, yeasts and other saprotrophs, was associated with low SOC stocks (F 1,68 = 17.41, p < 0.001, r 2 = 0.20; Figure 5d), whereas root-associated ascomycetes were more abundant on plots with high SOC stocks (F 1,68 = 6.74, p = 0.012, r 2 = 0.08; Figure 5d; Table S7).
When the individual and interactive impacts of grazing and tree vicinity were included in linear models of SOC stocks, fungal orders and functions explained SOC stocks even better (Table S8). For most taxa and guilds, accounting for the slightly higher SOC stocks in the vicinity of mountain birches (Figure 5a), that is, the hummocks formed around trees, the explanatory power of the model increased (Table S8). After including the effect of tree vicinity, also the abundance of root-associated basidiomycetes was negatively related with SOC stocks (Table S8).

| DISCUSS ION
We compared soil fungal communities in two contrasting longterm reindeer grazing regimes in an alpine treeline ecotone. In one regime, ~70 years of year-round grazing had hindered mountain birch establishment and led to an open landscape dominated by low stature evergreen dwarf shrubs, while the other regime with only wintertime grazing was characterized by a higher density of mountain birch (Oksanen et al., 1995). We found that the YGR was characterised by higher soil C:N ratio and that the shift in

| Linkages between fungal community composition and soil carbon storage
In line with our first hypothesis and some earlier studies linking reindeer grazing to altered fungal community composition and growth (Santalahti et al., 2018;Vowles et al., 2018), we found that root-associated ascomycetes, including ErM fungi, were most abundant away from mountain birches in the grazing regime exposed to year-round reindeer grazing and tightly linked to the abundance of ericaceous dwarf shrubs. In contrast, the abundance of root-associated basidiomycetes, mostly forming EcM, was linked to the relative abundance of Betula DNA in the soil.
Moreover, in support of our second hypothesis and observations from boreal forests (Clemmensen et al., 2013(Clemmensen et al., , 2015, we show that the abundance of root-associated ascomycetes correlated positively with soil carbon, whereas root-associated basidiomycetes correlated negatively with soil carbon when tree vicinityand thus hummock formation around tree trunks-was accounted for. However, in contrast to our third hypothesis and earlier observations from the treeline , the negative correlation of EcM fungi with soil carbon stocks could not be linked to EcM species that possess cord-forming mycelial structures. Instead, all the dominant EcM fungi in the study area, Russula claroflava, Russula decolorans, Russula paludosa, Pseudotomentella tristis (Thelephorales), Elaphomyces muricatus (Ascomyceta) and species from the order Cantharellales, were characterized by simple mycelia, whereas cord-forming EcM fungi contributed less than 1% of fungal communities. We suggest this pattern to be linked to the high sensitivity of cord-forming EcM, such as Cortinarius, to disturbances and climatic constraints (Vowles et al., 2018;Clemmensen et al., 2021). Thus, even the smaller disturbances induced by winter time grazing-the digging of soil and vegetation and the loss of an insulating snow layer in the grazing craters-could have benefitted EcM with shortdistance mycelia at the same time as the higher nutrient availability (i.e. low C:N ratio) reduced the demand for traits adapted to the break-down and long-distance transport of organically bound nutrients (Bödeker et al., 2014;Kyaschenko et al., 2017;Clemmensen et al., 2021).
Furthermore, in line with our hypotheses, we found that freeliving saprotrophic fungi were more abundant in the WGR and that their abundance correlated negatively with soil carbon storage.
These patterns were driven by the moulds Mortierellales, Mucorales, Eurotiales (Penicillium sp.) and the yeasts Leucosporidiales and Saccharomycetales. These moulds and yeasts are known for their opportunistic nature and dependence on easily available carbon as their own capacity to decompose complex OM is rather limited (Fernandez & Kennedy, 2018;Lindahl et al., 2010). Thus, the higher abundance of moulds and yeasts in the WGR is likely linked to the increased availability of easily decomposable leaf, root and mycelial litter (McLaren et al., 2017;Parker et al., 2018;Vankoughnett & Grogan, 2016). This is supported by the negative correlations between the abundance of saprotrophic fungi and DNA from plant groups with recalcitrant litters (i.e. bryophytes and the ericaceous species). As saprotroph abundance further correlated positively with EcM abundance, we suggest that they are particularly linked to mycelial biomass turnover of short-distance EcM fungi (Brabcová et al., 2016;Fernandez & Kennedy, 2018;Lindahl et al., 2010). In contrast to the ErM mycelia with melanised cell walls (Read & Perez-Moreno, 2003), EcM necromass is less recalcitrant and thus, easier to decompose by moulds and yeasts (Brabcová et al., 2016;Lindahl et al., 2010) and, in later stages of decomposition, also by EcM fungi (Fernandez & Kennedy, 2018). This activity of moulds, yeasts and high levels of allelopathic secondary compounds, which retard litter decomposition and nutrient mineralization rates (Tybirk et al., 2000). The recalcitrant necromass of root-associated ascomycetes may further push tundra heath ecosystems towards lower OM decomposition rates , not only because of the melanised cell walls of the ErM fungi (Fernandez & Kennedy, 2018) but also because the fungal necromass can form recalcitrant complexes with plant root derived condensed tannins, contributing to nutrients becoming locked up in soil OM (Adamczyk et al., 2019). Ultimately, these feedbacks could contribute to higher rates of soil carbon sequestration  as also indicated by the positive correlation between the abundance of root-associated ascomycetes and soil carbon stocks.

| Patterns in fungal diversity and abundance and the counterintuitive pattern in lichen abundance
The mountain birch dominated vegetation in the WGR supported higher species richness and effective diversity of the soil fungal community, whereas all diversity indices were indicatively lowest in the heath plots away from the trees in the YGR. These differences are likely linked to the high dominance of P. ericae and Archaeorhizomycetes in the soils away from the trees in the YGR.
Yet, the diversity pattern could also indicate a positive response of fungal diversity to an intermediate level of disturbance prevalent in the WGR. In contrast, fungal abundance (as inferred from ITS2 copies) was highest away from the trees in the YGR. Also this pattern is likely linked to the dominant taxa, P. ericae and Archaeorhizomycetes: the combination of the mycelias' high recalcitrancy (Fernandez & Kennedy, 2018) and limited decomposer capacity (Kyaschenko et al., 2017) could explain why higher fungal abundance did not contribute to lower soil carbon stocks (in contrast to Ylänne et al., 2020).
We found the abundance of DNA from lichenized fungi in the organic soil to be higher in the YGR compared with the WGR. This observation was contrary to our expectations, as both satellite images (Cohen et al., 2013;Forbes & Kumpula, 2009) and visual documentation indicate higher cover of particularly reindeer lichens (Cladonia spp.) in the WGR than in the YGR (32.1 ± 7.4% in WGR vs. 12.1 ± 3.9% in YGR). As also the lichenized DNA consisted mainly of Cladonia species, an important winter forage for reindeer, the results indicate a discrepancy between above-and belowground lichen abundance. This difference could derive from biases in the DNAbased data, such as differences in DNA extraction efficiency or variation in ITS2 copy numbers among species. Yet, given the known vulnerability of lichens to grazing and trampling (Bernes et al., 2015), we suggest that the more intense trampling in the YGR potentially resulted in a higher degree of fragmentation of the above-ground lichen communities and led to lichen necromass becoming incorporated into the litter and soil layers to a greater degree. The differences between the grazing intensities could further be accentuated by the impact of the environment on the decay rate of lichen necromass (Lang et al., 2009). As the abundance of lichenized fungi was particularly high under the trees in the YGR, where the OM rich hummocks underneath the trees were characterized by mosses from the families Polytrichum and Dicranum, these could contribute to a more acidic environment, where lichen fractions embedded in the moss layer remain, at least to a certain extent, protected from decay.

| Grazing impacts at the landscape level
When interpreting the results of this study at the landscape level, it is important to note that the sampling design largely overestimates the share of mountain birch habitats in the YGR and that the sampling approach did not include fungal communities in the deeper soil layers. Together, these factors could contribute to real differences between the grazing regimes being even larger than seen in this study. Due to differences in rooting depth of these vegetation types, also carbon stocks in deeper roots and soil may differ between the grazing regimes. Yet, the pattern in soil carbon across the landscape is unlikely to merely reflect the pattern in above-ground biomass and carbon uptake because trees and associated microbial communities also induce decomposition of old and more recalcitrant carbon (Hartley et al., 2012;Parker et al., 2020) resulting in an inverse relationship between the above-and belowground carbon stocks in mountain birch forests and tundra heaths (Hartley et al., 2012;Parker et al., 2015;Clemmensen et al., 2021). Our results highlight that such differences between these vegetation types may be induced by long-term reindeer summer-time browsing of deciduous trees. Yet, it remains unknown how common it is that grazers drive such vegetation transitions and how generalizable the belowground consequences are. In general, grazer impacts on ecosystems depend on the density, frequency and seasonality of the animal activity.
When occurring in high numbers in the snow-free period, grazers may limit deciduous shrub and tree establishment and growth (Biuw et al., 2014;Bråthen et al., 2017;Post & Pedersen, 2008 to certain extent, protect these from grazing (Christie et al., 2015).

| Increasing importance of grazing in the rapidly warming Arctic
It is well established that vertebrate herbivores, when occurring in high numbers and at a particular time, may exert negative effects on the abundance of deciduous shrubs and trees (Christie et al., 2015;Olofsson & Post, 2018). As vertebrate herbivory mainly prevents the establishment of new trees (Biuw et al., 2014;Kumpula et al., 2011), it only induces vegetation shifts from forest to open landscapes when accompanied by other disturbances that induce tree mortality, such as outbreaks of invertebrate herbivores (Biuw et al., 2014) or pathogens (Olofsson et al., 2011). However, ample evidence from across the arctic tundra has shown that vertebrate herbivores prevent climate-induced advancement of deciduous shrubs (Kaarlejärvi et al., 2015;Post & Pedersen, 2008;Zamin & Grogan, 2013) and keep tree-lines lower than their altitudinal limits (Cairns & Moen, 2004;Oksanen et al., 1995;Van Bogaert et al., 2011). In this context, it is becoming clear that vertebrate herbivory alone may be decisive for vegetation trajectories in a warming tundra (Bråthen et al., 2017;Christie et al., 2015;Olofsson & Post, 2018). This highlights the need to incorporate current grazing ranges for improved predictions of tundra vegetation development (Yu et al., 2017). Yet, it also emphasizes the need to predict the consequences of altered grazing patterns induced by societal and environmental changes, such as the dramatic declines in wild reindeer herds (Russell et al., 2018).
The northward and altitudinal expansion of shrubs and trees is occurring across Arctic and subarctic regions (Myers-Smith et al., 2019;Tape et al., 2012;Terskaia et al., 2020) and is expected to continue (Bjorkman et al., 2020). Although most observations of shrub expansion derive from riparian areas and tundra meadows (Bjorkman et al., 2020;Bråthen et al., 2017;Post & Pedersen, 2008), areas dominated by evergreen dwarf shrubs have also been reported to be encroached by tall deciduous shrubs in recent years (Terskaia et al., 2020). Such expansions are notable in scale as tundra heaths are the most abundant vegetation type in the Arctic, covering 13.6% of the vegetated land surface (classified as erect dwarf shrub tundra in the Circumpolar Arctic Vegetation map; CAVM Team, 2003). By combating the encroachment of deciduous shrubs and trees in these areas, both large and smaller herbivores could be decisive for shrub and tree ranges in a warmer climate. Here we suggest that when controlling the ranges of mountain birch forests and tundra heaths, herbivory has the potential to alter belowground fungal communities and soil carbon sequestration at the treeline ecotone. As these two vegetation types associate with different above-and belowground feedbacks linked to the dominant vegetation (Wardle et al., 2004) and to the soil fungal community (Read et al., 2004;Clemmensen et al., 2021), herbivory could control a tipping-point between the two vegetation trajectories and be decisive for the belowground consequences of climate change.

| CON CLUS IONS
Here we showed that ~70 years of different reindeer grazing regimes at the treeline ecotone shifted vegetation, soil fungal communities and soil carbon sequestration in disparate directions. We found that mountain birch woodlands experiencing only winter grazing hosted a fungal community characterized by short-distance EcM fungi and free-living saprotrophs and were linked to high soil fertility. In contrast, year-round grazing by reindeer promoted a more open tundra heath vegetation associated with a fungal community dominated by a few root-associated ascomycete species, whose abundance correlated with higher organic soil carbon stocks. These findings show that when controlling the abundance of deciduous shrubs and trees, vertebrate herbivores may have far-reaching consequences on key below-ground processes and, ultimately, on soil carbon sequestration. This supports earlier suggestions about the potential to use reindeer management (Bråthen et al., 2017) and even reintroductions of large grazers (Olofsson & Post, 2018) to reduce some of the unwanted effects of a warmer climate, such as shrub encroachment and the associated loss of soil carbon. Our study also highlights that past, current and future grazing patterns are important to include when predicting the belowground consequences of the ongoing climatic and environmental changes across the Arctic tundra.

ACK N OWLED G EM ENTS
The research presented in this paper is a contribution to the stra-

DATA AVA I L A B I L I T Y S TAT E M E N T
All data that support the findings of this study is openly available.