Materials and methods
Study design and sampling
We conducted sampling in September/October 2017 at nine different high-elevation forest line sites in south-central Norway (Fig. 1). Our selection of sites was based on a forest line model for Norway (Bryn & Potthoff, 2018). Based on a wider set of a priori candidate sites, we selected nine sites with a clear elevation gradient, a distinct shift from mountain birch forests to alpine vegetation and located at least 15 km from each other (Fig. 1a). At each site, we collected soil samples from plots located every 20 m along a 200 m transect, stretching from 100 m distance below the mountain birch forest line to 100 m distance above (i.e. 11 plots per transect and 99 plots all together; Fig. 1c). At each plot, five soil cores were collected in each orientation (i.e. North, East, South, West) and the centre within a circle with 1.5 m radius using a 3.8 cm in diameter soil corer, and pooled to one representative soil sample, excluding the mineral soil layer. Most often, there was no clear division between the litter and organic layers and we therefore analysed them together. We removed living green plants and mosses. Soil samples were kept at -20°C during fieldwork and later stored at -80 °C before further processing. Within each plot, we recorded the understory vegetation to species level for lichens, mosses and higher plants and categorized the abundance of each species into three levels: rare (1), common (2) and dominant (3). Based on this information, we calculated the proportion of each plant species per plot, as well as the overall proportion of the following plant groups: ErM plants, EcM plants, AM plants, lichens, mosses andPyrola species, the latter making arbutoid mycorrhiza.
We obtained the following site-specific environmental variables from a published data (Horvath et al., 2019): aspect, annual precipitation, mean temperature, bedrock and slope. The three types of bedrock provided in this dataset (“nutrient-poor bedrock”, “nutrient-average bedrock”, “nutrient-rich bedrock”) were ranked on a 1-3 scale. All environmental and climatic variables were then zero skewness transformed and standardized.
Soil analyses
Soil samples were thawed and handled further in a laminar flow hood, where plants and coarse (>2mm) pieces of wood and roots were removed. Approximately 70 g of soil sample was freeze-dried (Labconco corporation, Kansas City, MO, USA) in falcon tubes for 36 h and then pulverised using a FastPrep-24 beadbeater (M.P. Biomedicals, CA, USA) in two rounds of 20 sec at 25 MHz and subsampled for the different analyses. For soil pH measurements, we diluted 0.5 g of freeze-dried soil in 5 mL dH2O for one hour, and measured pH using a LAQUA-TWIN-11 pH Meter (Horiba Scientific, Kyoto, Japan) following the manufacturers protocol. Soil C and N concentration was determined, using 0.5 g of freeze-dried soil, by a flash elemental analyzer (Thermo Finnigan Flash EA 1112, ThermoFisher Scientific, Waltham, USA). Soil P concentration was determined by a segmented flow analyzer (SEAL AA3 HR AutoAnalyse, SEAL Analytical Ltd, Southampton, UK). We used approximately 200 mg of freeze-dried soil for measuring free and total soil ergosterol concentrations (mg g-1DW) using a similar protocol as in (Ransedokken, Asplund, Ohlson, & Nybakken, 2019).
Molecular methods
We extracted DNA using a CTAB-Chloroform DNA extraction protocol followed by a column based DNA purification using the E.Z.N.A.® Soil DNA Kit following the manufacturers protocol (Omega Bio-tek, Norcross, USA). Technical replicates and extraction negatives (negative controls) were introduced during DNA extraction, while mock communities (positive controls) were introduced during the PCR step. We PCR amplified both the rDNA ITS2 and 18S regions. The primers gITS7 (forward) and ITS4 (reverse) (Ihrmark et al., 2012) were used for amplifying the ITS2 region, and TAReuk454FWD1 (forward) and TAReukREV3 (reverse) (Stoeck et al., 2010) for 18S. All primers were tagged with unique molecular identifiers (MID). Each PCR reaction consisted of 1 µl DNA template and 24 µl master mix: 15.7 µl dH2O, 2.5 µl Gold Buffer, 2.5 µl Gold MgCl2, 1 µl 20mg/ml BSA, 0.2 µl dNTPs, 0.1 µl AmpliTaq Gold, 1.5 µl 10 µM forward primer and 1.5 µl 10 µM reverse primer. We run PCR reactions for ITS2 with initial denaturation at 95°C for 5 min, followed by 32 cycles of denaturation at 95°C for 30 sec, primer annealing at 55°C for 30 sec and elongation at 72°C for 1 min. We added a final elongation step at 72°C for 7 min, before cooling down to 4°C. The protocol for 18S was slightly different: PCR reactions were run with initial denaturation at 98°C for 7 min, followed by 32 cycles of denaturation at 98°C for 30 sec, primer annealing at 53°C for 30 sec and elongation at 72°C for 45 sec. A final elongation step was included at 72°C for 10 min, before cooling down to 4°C. We controlled each PCR product for positive amplification with gel electrophoresis using a 2% agarose gel, before individual clean-up and purification of the amplicons with ZR-96 DNA Clean & Concentrator-5 kit (Zymo Research, California, USA). DNA concentrations for each sample were measured with the Qubit 2 fluorometer dsDNA BR Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA) and pooled to equimolar concentration into four pools. Each pool was cleaned and concentrated with DNA Clean & Concentrator-5 (Zymo Research, California, USA). The four libraries (two for each targeted DNA region) were sequenced by Fasteris SA (Switzerland) in two runs using Illumina MiSeq with a 250 bp paired-end (PE) with V3 chemistry. A ligation protocol, specifically designed to minimize tag-jumping, was used to ligate the amplicons with the Miseq flow-cell adapters.
Bioinformatics
We performed all bioinformatics analyses on the Abel high-performance computer cluster at the University of Oslo. The raw PE reads were demultiplexed separately with simultaneous tags and primers removal using CUTADAPT (Martin, 2011). No miss-match with primer and MID tags were allowed. Further processing of the data was performed with the DADA2 pipeline (Callahan et al., 2016) using the statistical environmentR version 3.6.0 (R Core Team, 2014). We generated sequence quality profiles and used this information to decide parameters for filtering low quality reads. A truncate length of 200 bp was set for the 18S sequence reads (hereafter called reads). Due to natural length variability of the ITS2 marker, no truncation was imposed. In both datasets, maximum expected error was set to 2.5. We then independently corrected forward and reverse reads using the DADA2 machine-learning algorithm that estimates correction parameters from the data itself. The default settings were used for 18S, while for ITS2 the band size was set to 32 due to the higher INDEL rate of this marker. The corrected forward and reverse reads were merged. We used a minimum overlap of 50 and 8 nucleotides when merging ITS2 and 18S reads, respectively. Chimeras were checked and removed with a denovo approach using the DADA2 bimera algorithm with default settings. We then constructed an amplicon sequence variant (ASV) table with the chimera-free reads. For the ITS2 data, we used ITSx (Bengtsson-Palme et al., 2013) for extracting the ITS2 region and filtering of non-fungal reads. Due to widespread intraspecific variation in ITS2, an additional clustering step with 97% similarity was performed using VSEARCH (Rognes, Flouri, Nichols, Quince, & Mahe, 2016) and singleton sequences were discarded. Finally, to adjust for over-splitting of Operational Taxonomic Units (OTUs), we performed post clustering curation using the LULU algorithm (Froslev et al., 2017) for both datasets.
For taxonomic annotation, the ITS2 dataset was query searched using VSEARCH global against UNITE v8.0 (Nilsson, Larsson, et al., 2019), whereas the 18S dataset was annotated using the Eukref/PR2 database (del Campo et al., 2018). We used FunGuild (Nguyen et al., 2016) for functional annotation of the ITS2 dataset, where we annotated the fungi into the following functional groups: Saprotrophs, ectomycorrhizal fungi (EcM), lichens, yeasts, pathotrophs, root-associated ascomycetes (including DSEs and ErM fungi) and other ascomycetes. Notably, these groups should be regarded tentative (i.e. as hypotheses), since we still lack solid evidence for most fungi’s ecologies and functions.
In the ITS2 dataset, we removed one sample due to low read number (286 reads). Number of reads per sample was rarefied to the lowest sample read number (6317 reads). The final ITS2 dataset consisted of 98 samples and 3090 OTUs. In the 18S dataset, we removed six samples due to low read number (< 9). Plant OTUs were removed, and number of reads per sample was then rarefied to lowest sample read number (2486 reads). The final 18S dataset consisted of 93 samples and 4595 OTUs.
Statistics
We performed all statistical analyses in R (R Core Team, 2014).ggplot2 was used for graphical plotting (Wickham & Chang, 2016); the vegan package for community composition analyses, i.e., rarefaction, OTU count data transformation, ordinations, variation partitioning, environmental correlations with ordination analyses (Oksanen et al., 2007) and the nlme (Pinheiro et al., 2017) as well as MuMIn (Barton, 2019) packages for modelling correlations.
To investigate structure in community composition, we constructed ordination diagrams using global non-metric multidimensional scaling (GNMDS) with the “metaMDS” function on the rarefied OTU matrix, using settings as recommended by Liu et al. (Liu et al., 2008). GNMDS and detrended correspondence analysis (DCA), using default settings, were run in parallel to assess the correlation between GNMDS and DCA axes. As shown by Son and Halvorsen (van Son & Halvorsen, 2014), a good correspondence between these methods can be interpreted as robust ecological gradients. We fitted the environmental variables to the ordinations with the “envfit” function, and plot isolines were fitted with the “ordisurf” function. In order to investigate the proportion of variation explained by different environmental variables, variation partitioning of the ITS2 and 18S data was performed using canonical correspondence analysis (CCA) function in the vegan package. We ran a forward model selection until no more significant variables could be added. Correlations between soil C and ergosterol and DSEs were modelled as linear mixed models with sites as a random factor using the “lme” function, and r squared values extracted from the models using the “r.squaredGLMM” function.