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.