Abrupt changes in the composition and function of fungal communities along an environmental gradient in the high Arctic

Fungi play a key role in soil–plant interactions, nutrient cycling and carbon flow and are essential for the functioning of arctic terrestrial ecosystems. Some studies have shown that the composition of fungal communities is highly sensitive to variations in environmental conditions, but little is known about how the conditions control the role of fungal communities (i.e., their ecosystem function). We used DNA metabarcoding to compare taxonomic and functional composition of fungal communities along a gradient of environmental severity in Northeast Greenland. We analysed soil samples from fell fields, heaths and snowbeds, three habitats with very contrasting abiotic conditions. We also assessed within‐habitat differences by comparing three widespread microhabitats (patches with high cover of Dryas, Salix, or bare soil). The data suggest that, along the sampled mesotopographic gradient, the greatest differences in both fungal richness and community composition are observed amongst habitats, while the effect of microhabitat is weaker, although still significant. Furthermore, we found that richness and community composition of fungi are shaped primarily by abiotic factors and to a lesser, though still significant extent, by floristic composition. Along this mesotopographic gradient, environmental severity is strongly correlated with richness in all fungal functional groups: positively in saprotrophic, pathogenic and lichenised fungi, and negatively in ectomycorrhizal and root endophytic fungi. Our results suggest complex interactions amongst functional groups, possibly due to nutrient limitation or competitive exclusion, with potential implications on soil carbon stocks. These findings are important in the light of the environmental changes predicted for the Arctic.

Fungi are essential for the functioning of arctic ecosystems and control much of the carbon (C) balance (Dahlberg & B€ ultmann, 2013).
Root-associated fungi help vascular plants to obtain N, which is highly limiting in arctic soils (Hobbie & Hobbie, 2006). Saprotrophic fungi are active primary or secondary decomposers and strongly influence C and nutrient cycles (Christiansen et al., 2016), while pathogenic fungi can cause severe diseases in plants and can greatly decrease plant biomass (Olofsson, Ericson, Torp, Stark, & Baxter, 2011). Many arctic fungi form lichens, that are important primary producers (Longton, 1988;Webber, 1974), while also contributing to N fixation (Crittenden & Kershaw, 1978), and are important food sources for several herbivores. However, the abiotic and biotic factors that determine the occurrence of these functional groups at landscape (e.g., habitat) and microhabitat scale and the link between composition and function are still poorly understood, especially for nonlichenised fungi (Dahlberg & B€ ultmann, 2013). Several studies have focused on the abiotic and biotic factors that influence the distribution of specific functional groups in the Arctic. For example, the richness and composition of ectomycorrhizal (ECM) fungal communities covaried with landscape scale differences in abiotic conditions (due to changes in snow depth, soil nutrient content or pH) and biotic conditions (due to changes in plant richness, floristic composition or vegetation cover) (Fujimura & Egger, 2012;Morgado et al., 2016;Semenova et al., 2016;Timling et al., 2014). In addition, the biotic interactions occurring at microhabitat scale (e.g., between hostdependent fungi and Dryas or Salix shrubs, Fujimura and Egger (2012)) are also expected to shape the composition of some groups, such as ECM fungi. The influence of these biotic interactions on the composition of ECM fungal communities at microhabitat scale remains, however, unclear. Some authors report a strong specificity with host plant (Fujimura & Egger, 2012), whereas many others found no clear specificity (Blaud, Phoenix, & Osborn, 2015;Botnen et al., 2014;Mundra et al., 2014;Ryberg, Larsson, & Molau, 2009;Walker et al., 2011). Furthermore, the composition of fungal communities can also be shaped by the interaction amongst functional groups, as described for ECM and saprotrophic fungi, because they likely compete for nutrients (Leake, Donnelly & Boddy, 2002;Nordin, Schmidt, & Shaver, 2004;Wallenstein et al., 2007).
Only very few studies, however, have analysed fungal communities accounting for several co-occurring functional groups along environmental gradients in the Arctic (but see Timling et al., 2014;Geml et al., 2015 andSemenova et al., 2016). These studies indicated that the composition of the main functional groups of fungi in the Arctic (i.e., saprotrophic, lichenised, plant pathogenic, root endophytic and ECM fungi) can be structured at different spatial scales, from bioclimatic zone to habitat or even microhabitat scales in response to environmental variability.
In this study, we analyse the impact of habitat and microhabitat on the taxonomic and functional composition of fungal communities in Northeast Greenland. We investigated soil fungal communities in three contrasting habitats (e.g., the snowbed, heath and fell field) that are widespread across mesotopographic gradients in the Arctic and that are associated with an increase in environmental severity from the snowbed to the fell field (Walker, 2000). In each habitat, we analysed the soil fungal communities occurring in three widespread microhabitat types: (i) patches dominated by the dwarf shrub Dryas octopetala L. 9 intermedia Vahl., (ii) patches dominated by the dwarf shrub Salix arctica Pall. and (iii) patches without a dominating shrub, mostly bare ground, or with a very thin cryptogamic layer (bare ground patches hereafter). Both Dryas and Salix are dominant ECM hosts in the high Arctic (Newsham, Upson, & Read, 2009;. We hypothesized that: (i) the differences in the composition of fungal communities would respond primarily to habitat type, (ii) the composition of fungal communities would be sensitive to microhabitat type within each habitat, (iii) the proportional richness of the functional groups (i.e., saprotrophic, lichenised, plant pathogenic, root endophytic and ECM fungi) would covary with abiotic conditions, and (iv) groups with intimate interactions with vascular plants (ECM, root endophytic and plant pathogenic fungi) would dominate in the Dryas and Salix patches. To our knowledge, this study is the first to use high-throughput molecular methods to explore soil fungal communities in Greenland and one of the first conducted in the Arctic that studies several co-occurring functional groups of fungi along a natural environmental gradient.

| Study area and data collection
The fieldwork was conducted in July 2011 at Zackenberg Valley (Northeast Greenland; 74°28'N, 21°33'W). The climate in this area is strongly affected by the coastal ice belt, which provides a continental climate, with very cold winters, little precipitation and sunny summers (Meltofte & Rasch, 2008). Mean monthly temperatures in winter are below À20°C, and the wind is mostly from the north. The mean monthly temperatures in the short snow-free summer usually vary between 3 and 7°C in July and August. The mean annual precipitation is 250 mm (Hansen et al., 2008). Snow covers the ground for a minimum of 6-7 months; the growing season starts in late May in early snow-free areas, and the snow cover is extensive until early summer in snowbeds (Hinkler, 2005;Meltofte & Rasch, 2008). Automatic snow measurements conducted from 2005 to 2015 indicate that strong winds during the cold season tend to accumulate the snow at low elevations and leeward sides, whereas the snow almost disappears from the exposed ridges (Elberling et al., 2008;Hinkler, 2005;Walker, 2000). Snow depth thus generally decreases from low to high elevation. GRAU ET AL.

| 4799
We chose three different habitats along a topographic gradient on the southwestern slope of Mount Aucellabjerg associated with an increase in environmental severity: (i) a snowbed at ca. 40 m a.s.l., (ii) a heath at ca. 200 m a.s.l and (iii) a fell field at ca. 425 m a.s.l. Snow depth in the snowbeds is 1-1.5 m in winter and ranges from 0 to 30 cm in the fell field (http://zackenberg.dk/publications/annua l-reports/, accession date: July 2016). Snow depth was not automatically measured in the heath, but direct measurements indicated that snow depth usually reached 60 cm in winter (Elberling et al., 2008).
These measurements therefore confirm the large, moderate and very shallow winter snow cover in the snowbed, heath and fell field, respectively, as expected from the topography. At landscape scale throughout the Arctic, topography drives drainage and soil moisture, which, in turn, greatly influence soil development, the composition of tundra communities and ecosystem processes (Walker, 2000). For example, Elberling et al. (2008) described the abiotic conditions at Zackenberg region and reported that (i) the depth of the active soil layer increased from the snowbed to the fell field, (ii) soil-water saturation at the soil surface decreased from 65% to 90% in the snowbed to 40%-60% in the heath and to <40% in the fell field, (iii) soil nutrient content decreased from the snowbed to the fell field.
These differences in abiotic conditions across habitats produce a gradient of increasing environmental severity from the snowbed to the fell field (referred hereafter as the "main gradient") and a corresponding decrease in plant cover (Table S1).
We selected three replicate 10 9 10 m plots in each habitat at similar elevation and aspect, separated by a few hundred metres ( Fig. S1). In each plot, we selected patches larger than 25 9 25 cm where Dryas, Salix, or bare ground clearly dominated (with >80% cover). Four soil samples were collected for each microhabitat in each plot with a 3.5-cm corer to a depth of 5 cm. The coarse litter layer was removed in the patches of Dryas or Salix, and the soil cores were centred on the base of the roots to collect soil from the rhizosphere. The four soil samples were mixed and pooled into one composite sample (n = 3 habitat types 9 3 microhabitats 9 3 plots = 27 composite samples in total). The samples were immediately frozen for transport to the laboratory and were divided into two subsamples: one was preserved in tubes with cetyltrimethylammonium bromide (CTAB) for subsequent molecular analyses, and the other was used to measure pH and the contents of organic C, total N, and available phosphorus (P). Organic carbon was quantified by acid dichromate oxidation, following the modified Mebius' method (Nelson & Sommers, 1996), with the added modification that oxidation was performed at 155°C, as recommended by Soon and Abboud (1991). Total nitrogen was analysed in a ThermoQuest CHN analyser. Available phosphorus was analysed following (Olsen & Sommers, 1982). Soil pH was measured with a glass electrode, in 10:25 soil-water suspensions (weight volume), except in a few soil samples very rich in organic matter, in which a proportion 10:50 was applied.
We also assessed the floristic composition where the soil samples were collected. See Grau, Ninot, P erez-Haase, and Callaghan (2014) for further details on the floristic composition and vegetation patterns in the study plots.

| Molecular analyses and bioinformatics
DNA was independently extracted twice from 12 ml of suspensions of each sample containing CTAB and at least 5 g of dried soil using a combination of a CTAB-based method and the method used in the Macherey-Nagel (Gmbh & Co., D€ uren, Germany) soil kit. The soil suspension in CTAB was incubated at 55°C for 15 min and then centrifuged for 5 min at 12,000 g to pellet the soil debris. The supernatant was thoroughly mixed with 24:1 chloroform:isoamyl alcohol and centrifuged at 13,000 rpm for 1 min. The upper aqueous phase containing the DNA was mixed with an equal volume of binding buffer (Macherey-Nagel) and filtered through a DNA-binding silica membrane. DNA binding and purification followed the Macherey-Nagel protocol. This method produced highly concentrated and pure DNA. The two extractions were pooled prior to PCR. The ITS2 region (ca. 250 bp) of the nuclear ribosomal rDNA repeat was PCR amplified using primers fITS7 (Ihrmark et al., 2012) and ITS4 (White, Bruns, Lee, & Taylor, 1990), as described in detail in Geml et al. (2014). The ITS4 primer was labelled with sample-specific Multiplex Identification DNA tags (MIDs). The amplicon library was sequenced at Naturalis Biodiversity Center (Leiden, The Netherlands) using an Ion 318TM Chip and an Ion Torrent Personal Genome Machine (PGM; Life Technologies, Guilford, USA).
The raw data (1625417 sequence reads) were initially filtered using the online platform Galaxy (https://main.g2.bx.psu.edu/root), with which the sequences were sorted by sample, and adapters and MIDs were removed. Primer sequences were removed, and poor-quality ends were trimmed off based on a 0.02 error probability limit by Geneious Pro 5.6.1 (BioMatters, Auckland, New Zealand). The sequences were then filtered using MOTHUR v. 1.32.1 (Schloss et al., 2009) using the following settings: no ambiguous bases, homopolymers no longer than 10 nucleotides, and length range from 150 to 400 bp, which provided 1195 838 quality-filtered sequences with an average read length of 253.3 AE 47.4 (mean AE SD) bp. Identical sequences were collapsed into unique sequence types while preserving corresponding read counts, and singletons (568743; 47.56%) were removed with USEARCH v.8.0 (Edgar, 2010). The remaining sequences were grouped into 1911 operational taxonomic units (OTUs) at 97% sequence similarity using USEARCH while simultaneously removing putative chimeric sequences (6851; 0.6%). We assigned sequences to taxonomic groups based on pairwise similarity searches against the curated UNITE fungal ITS sequence database (accessed in 2015) containing identified fungal sequences with assignments to Species Hypothesis groups (Kõljalg, Nilsson, Abarenkov, & Al, 2013). After excluding OTUs with < 80% similarity or < 150-bp pairwise alignment length to a fungal sequence, 1601 fungal OTUs were retained. The representative sequences of fungal OTUs have been submitted to GenBank. These OTUs were then assigned to functional groups based on the identities of the most similar sequences (above 90% sequence similarity to avoid uncertainty), coupled with the source information of the reference sequence, when available, or based on published ecological information of the taxon in question (e.g., Nguyen et al., 2016; Tedersoo & Smith, 2013). The data set containing only fungal OTUs comprised a total of 610345 high-quality sequences with an average of 22605 AE 7400 reads per sample. Next-generation sequencing libraries are generally normalized by random subsampling to the size of the smallest library, but this practice often preserves only a small fraction of the useable data. This can increase type I (decreased specificity) and type II (decreased sensitivity) errors and was strongly discouraged by McMurdie and Holmes (2014). To avoid this problem, we opted to include all high-quality fungal sequences in the final analyses and we only used rarefied data for exploratory purposes and to check the robustness of our data by rarefying the sequence counts of all samples to the lowest library size (6281 sequences). We also calculated rarefied OTU accumulation curves ("vegan" R package; Oksanen et al., 2015) to explore the completeness of our sampling.
The sequencing depth of our samples was very high and rarefaction curves reached the saturation point in all samples (Fig. S2).

| Data analysis
We retained the OTUs with five or more sequences (1177 out of 1601 OTUs) to minimize false positives and contaminants with low abundances, as suggested by Lindahl et al. (2013). To avoid the uncertainty whether read abundance was a good indicator of OTU abundance in the samples (Amend, Seifert, & Bruns, 2010), we transformed the read-abundance matrix into a presence-absence matrix.
To investigate the similarity of the fungal communities amongst habitats and microhabitats, we ran nonmetric multidimensional scaling (NMDS) ordinations with the presence-absence matrix using the Jaccard dissimilarity index and with the read-abundance data using the Bray-Curtis index ("BiodiversityR" package; Kindt & Coe, 2005).
Because the stress value of the NMDS was slightly higher when read-abundance data were used (0.092 vs 0.072), we focused on the presence-absence data. For exploratory purposes and to check the robustness of the data, we also performed ordinations using rarefied data. We calculated the proportional richness of various functional groups (number of OTUs of a given functional group/total number of identified OTUs on a per-sample basis) to assess and compare their prevalence amongst the habitats and microhabitats and to account for the variation in the number of OTUs assigned to the functional groups amongst the samples (63% of all OTUs). This percentage is similar to or higher than that found in studies on the functional ecology of fungi (e.g., Semenova et al., 2015;Tedersoo et al., 2014) and we deem our conservative and reliable assignment appropriate to describe the general ecological patterns of the most common functional groups. We used the "envfit" R function to fit these proportions, the soil variables, fungal richness, and plant richness onto the NMDS ordination. We tested the effect of habitat and microhabitat by a nested PERMANOVA using the "adonis" and "nested.npmanova" functions; habitat was a fixed factor, microhabitat was a fixed factor nested within plot, and plot was a random, blocking factor (two-way split-plot design; Fig. S1). In addition, we used partial Mantel tests in PC-ORD v. 6.0 (Bruce & Grace, 2002) to differentiate the effects of floristic composition and abiotic variables on community structure. For the latter, mean and standard deviation values of the environmental variables were standardized so that all variables could have the same influence on the data. The data on floristic composition were obtained from Grau et al. (2014). We also ran an NMDS separately for each functional group to analyse the within-group variation in fungal composition amongst habitats. We did not analyse the effect of microhabitat on each functional group separately to simplify the analysis. We performed pairwise t tests with Bonferroni correction to identify differences in abiotic conditions, total fungal richness, and proportional richness of functional group amongst habitats and microhabitats. We assessed the occurrence of indicator OTUs in each habitat and microhabitat by a multilevel pattern analysis ("multipatt" function, "indicspecies" R package;  Table S2). The fell field had the most indicator OTUs.
The richness pattern described above was apparently driven by the phylum Ascomycota (617 OTUs; Fig. S3a), whereas the richness of Basidiomycota (410 OTUs) did not differ amongst habitats (Fig. S3b).
The proportional richness of Ascomycota also decreased significantly towards the snowbed, while that of Basidiomycota increased (Fig. S3c, d). The most common orders of Ascomycota and Basidiomycota are shown in Fig. S4. Both the total richness and the proportional richness of saprotrophic, lichenised and plant pathogenic fungi tended to decrease from the fell field to the snowbed (Figure 1b-d, Table S3).
The total richness of root-associated fungi (ECM fungi and root endophytes), however, did not differ significantly amongst habitats, although their proportional richness increased in the heath and in the snowbed (Figure 1e-f) due to the decrease in total fungal richness.
Overall, the proportional richness of ECM fungi and root endophytes was, therefore, negatively correlated with that of saprotrophic fungi, plant pathogenic fungi and lichenised fungi (Fig. S5). When the microhabitats were analysed separately within each habitat, the total fungal richness and the proportional richness of the functional groups in most cases did not differ much amongst Dryas, Salix, and bare ground (Figure 1g-l, Table S3). In parallel, the proportional richness of the functional groups in a given microhabitat was variable across habitats (e.g., Salix in the fell field vs Salix in the heath or in the snowbed). GRAU ET AL. | 4801

| Community composition at habitat and microhabitat scales
The fungal community composition of OTUs based on presence-absence data varied significantly across habitats and to a lesser extent by microhabitat type (Figure 3a). PERMANOVA analyses indicated that compositional turnover was significant amongst habitats (p < .005) and also amongst microhabitats within each habitat (p < .001). When read-abundance data were analysed by NMDS, we observed essentially the same patterns in response to habitat and microhabitat types; the patterns were also very similar whether we used rarefied or nonrarefied data (see Fig. S6a-c). Partial Mantel tests indicated that both abiotic variables and floristic composition correlated significantly with fungal community composition. However, when abiotic factors were controlled, the effect of floristic composition on fungal community composition was weaker (r = .314213, p = .000625) than the effect of abiotic factors (r = .740132, p = .000063) when floristic composition was accounted for (control matrix).

| Correlation between soil variables and compositional and functional variability
Total contents of N and organic C were highest in the snowbed, the content of available P increased progressively from the fell field to the snowbed, and the C:N and pH decreased towards the snowbed  Table 1) but not with plant or fungal richness. The within-group variation was explained mostly by habitat type and thus by the varying environmental conditions along the main gradient ( Figure 6). Information for the regression weights for all identified fungal families is also available in Table S4, but the weights are not discussed further in this report. Information on the Species Hypothesis for each OTU, the Accession no. in GenBank and their taxonomic and functional assignment can be found in Table S5.

| DISCUSSION
The data presented here show that both fungal richness and community composition are primarily driven by habitat-scale differences in abiotic conditions along the mesotopographic gradient studied, which supports our hypothesis 1. For example, one of the more striking patterns is that the severity of the environment is positively correlated with the total fungal richness. In other words, the most hostile environment (fell field) apparently hosts the highest number of fungal species. However, richness patterns of functional groups of fungi differ greatly along the gradient, as discussed below. In addition to habitat, the microhabitat type also has a significant, although less strong, effect on community composition, supporting hypothesis 2 as well (     (Geml, Kauff, Brochmann, & Taylor, 2010;Geml et al., 2012a,b), the richness and composition of the established species pool at any arctic location is primarily the result of coarse-scale environmental filters, such as regional edaphic and climatic conditions (Cox, Newsham, Bol, Dungait, & Robinson, 2016;Fierer et al., 2012;Goldmann et al., 2016;Tedersoo et al., 2014). It is noteworthy that fungal community in the fell field was not only the richest, but the compositionally most unique (Figures 2   and 3) and contained the most indicator species (Table S2) we also observed compositional differences amongst habitats within the functional groups ( Figure 6). This suggests that at a habitat scale, within-group functional variability (e.g., differences in nutrient acquisition or contrasting types of hyphal exploration) may influence species composition in relatively stable, climax communities by nichebased rather than stochastic processes that are more common in early successional stages (Blaalid et al., 2012;Jumpponen, 2003).

| Functional responses of fungal communities to environmental conditions at the habitat scale
The increase in environmental severity from the snowbed to the fell field was manifested by a decrease in snow depth, water saturation and nutrient content in the soil (Table S1, Figure 4). Soil temperature during the cold season is also expected to covary with snow depth (Brown, Hinkel, & Nelson, 2000;Woo, Mollinga, & Smith, 2007), with lower soil temperatures in areas with little snow accumulation (Jones, 1999), thereby contributing to the increase in environmental severity towards the fell field. This is also illustrated by the decrease in vegetation cover (87%, 85%, and 28% in the snowbed, heath and fell field, respectively; Table S1).
The decrease in fungal richness towards the snowbed may have been associated with the decrease in richness in saprotrophic, lichenised, and plant pathogenic fungi (Table S3) F I G U R E 4 Total N content, available P content, C:N, pH, and organic C content in each habitat (a-e) and microhabitat type (f-j; dark grey, Dryas; light grey, bare ground; white, Salix). The letters in a-e indicate significant differences in pairwise t tests with Bonferroni correction (significant, p < .05) been observed in habitats with thick snow cover . The changes in the relative importance of each functional group were closely associated with the environmental severity (Figures 1, 3b and 5). The proportional richness of root-associated fungi was higher in the snowbed and heath, likely due to more favourable abiotic conditions and increased abundance of their hosts (Table S1; Grau et al., 2014). High shrub cover generally corresponds to many roots available for mycorrhizal colonization, and mycorrhizal fungi favour the development of dominant shrubs such as Dryas and Salix .
The proportional richness of mycorrhizal and saprotrophic fungi had opposite trends along the main gradient. Mycorrhizal plants are strong competitors for N, as they can short-circuit the microbial mineralization step and compete against microbial decomposers for Nrich compounds (Nordin et al., 2004). Plants and their mycorrhizal symbionts can outcompete saprotrophic fungi and other microbial decomposers by decreasing the availability of N in the soil (Leake, 2002), which may account for the opposite patterns in mycorrhizal and saprotrophic fungi along the main gradient (Figures 1 and 3b, Table S3). The decrease in decomposers, due to competitive exclusion by mycorrhizal fungi, is expected to slow soil C respiration and increase soil C storage (Orwin, Kirschbaum, St John, & Dickie, 2011), as demonstrated recently by Averill and Hawkes (2016). This mechanism could, therefore, explain the higher accumulation of soil organic C in the snowbed and the heath (Figure 4). Moreover, mycorrhizal fungi store large amounts of C in their hyphae (Clemmensen et al., 2013), which may also contribute to the high accumulation of C in habitats with a high proportional richness of root-associated fungi.
In contrast, competitive exclusion by mycorrhizal fungi should be less pronounced in the fell field, because the proportional richness of mycorrhizal fungi was much lower. As a consequence, the higher proportional richness of saprotrophic fungi possibly explained the lower accumulation of organic C in the fell field (Figure 4), together with the very low input of dead biomass in this poorly vegetated habitat (Table S1). The increase in proportional richness of lichenised

Lichenised (%)
F I G U R E 5 Scatter plots between the response variables (total fungal richness and the proportional richness of saprotrophic, lichenised, plant pathogenic, root endophytic and ectomycorrhizal fungi) with abiotic variables (total N content, available P content, C:N, pH and organic C content). The fitted lines represent Loess smoothers. All regressions were statistically significant (p < .05) when analysed individually in a linear model (except for those indicated with n.s., not significant or m.s, marginally significant, p<.1) fungi towards the fell field (Figures 1 and 3b, Table S2) is likely due to the decrease in vegetation cover, because lichens are easily outcompeted by dense vascular vegetation but thrive in poorly vegetated habitats (Cornelissen et al., 2001;Joly, Jandt, & Klein, 2009).
The higher proportional richness of pathogenic fungi in the fell field is surprising at first glance given the low plant cover. However, most plants in this hostile habitat likely are closer to their physiological limits and may be more prone to pathogenic infections (Read, 2003).
The functional patterns were very similar to those reported by

| Composition and function of fungal communities at a microhabitat scale
The within-habitat differences in the composition of fungal communities were not explained by any of the abiotic variables measured, because the abiotic conditions did not vary consistently amongst the microhabitats (Figure 4 f-j). The presence of Dryas or Salix may hypothetically cause small microscale differences in abiotic variables that were not measured, such as soil temperature or water content, but such small differences alone would not likely be able to explain the differences in composition. The lack of clear abiotic differences may account for the lack of differences in the proportional richness of the functional groups amongst microhabitats.
We expected that ECM, root endophytic and plant pathogenic fungi would have a higher proportional richness in Dryas and in Salix patches and not in bare ground due to the biotic interaction expected to occur between these fungal groups and Dryas or Salix, temperatures, deeper snow cover in winter (Hinkler, 2005;Stendel et al., 2008) and earlier spring snow melt (Schmidt, Kristensen, Michelsen, & Bay, 2012). These climatic trends are consistent with the proliferation of shrubs in the Arctic (Sturm, Racine, & Tape, 2001;Wookey et al., 2009). The composition and function of fungal communities are expected to shift markedly with the changes in abiotic conditions and the expansion of shrubs (Dahlberg & B€ ultmann, 2013;Dumbrell, Nelson, Helgason, Dytham, & Fitter, 2010). This expectation is consistent with the clearly different composition and functions of the fungal communities amongst habitats with contrasting environmental conditions in our study. This is important for the understanding of potential changes in fungal composition and function (Peay, Bidartondo, & Arnold, 2010) under contrasting scenarios of environmental conditions.
If future abiotic conditions favour the expansion of habitats with high shrub cover towards poorly vegetated habitats, there could be marked changes in compositional and functional characteristics of fungal communities. Future warming-induced changes could also potentially decrease the extension of fell fields and cause local extinctions of the particularly unique set of fungi occurring in this habitat. We predict that mycorrhizal communities will continue to dominate in habitats where shrub cover is already very high, but the composition of fungal species may shift, especially if there are changes in the floristic composition of plant communities (Dahlberg & B€ ultmann, 2013). We hypothesize that the proportion of saprotrophic fungi will decrease if mycorrhizal communities are favoured, similarly to what we observed in the habitats where root-associated fungi prevailed, with potential impacts on nutrient cycling and C flow (Averill & Hawkes, 2016;Orwin et al., 2011;Wallenstein et al., 2007). Terricolous lichens are expected to face increasing competition from vascular plants, which may in turn have a negative impact on herbivores that feed on lichens (Joly et al., 2009).

DATA ACCESSIBILI TY
The list of OTUs identified in this study, their Accession no. in Gen-Bank and the taxonomic and functional assignments are available in the Supplementary Material (Table S5)