Here we present the data and R script from “Reindeer control over subarctic treeline alters soil fungal communities with potential consequences for soil carbon storage” by Henni Ylänne, Rieke L. Madsen, Carles Castaño, Daniel B. Metcalfe and Karina E. Clemmensen (Global Change Biology, 2021). In this study we reported the impacts of grazing regime and mountain birch vicinity on the abundance, diversity and community composition of the soil fungal community, and explored how the soil fungal communities relate to vegetation and, ultimately, to the observed variation in soil carbon stocks. The study was conducted at a treeline ecotone in northern Fennoscandia, where two adjacent 65-year-old grazing regimes with or without summer grazing by reindeer (Rangifer tarandus) could be compared. As these grazing regimes differ in the abundance of mountain birch (Betula pubescens ssp. czerepanovii), the own effects of mountain birch vicinity on the fungal communities were assessed by collecting soil samples directly underneath and > 3 m away from mountain birches.

This dataset includes data of soil properties, soil carbon and nitrogen stocks, vegetation community composition (as deterred from plant reads captured by the high-throughput sequencing of ITS2 region), fungal abundance (ITS2 copy numbers), fungal diversity and the taxonomic and functional composition of the soil fungal communities.

Study site, experimental design and soil sampling

The study was conducted at Jávrrešduoddarat, an alpine treeline ecotone at the border between Finland and Norway (68°45ʹN, 23°43ʹE; altitude 430–480 m a.s.l.). Since the late 1950’s, the Norwegian side of the Jávrrešduoddarat area is grazed by reindeer only during snow-covered winter months. In contrast, the Finnish side of the Jávrrešduoddarat area is grazed all year round. As reindeer exert less control over vegetation during winter months, when soil and vegetation are mostly protected by snow and ice, the year-round grazed side is more exposed to reindeer browsing, trampling and fertilization by excreta than the winter grazed side. This long-term grazing difference is reflected in higher lichen cover and greater abundance of mountain birch (Betula pubescens ssp. czerepanovii) in the winter grazing regime when compared to the year-round grazing regime.

Here we compared the open woodland in the year-round grazing regime (YGR) to the denser mountain birch woodlands in the winter grazing regime (WGR) by selecting three blocks (A–C), 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. In each block, six coordinates were marked on either side of the fence, 20 m from the fence with 50 meters between neighbouring points. From all coordinates, we collected soil samples directly underneath and > 3 m away from mountain birches, amounting to a total of 72 soil samples (3 blocks × 2 grazing regimes × 6 coordinates × 2 distances to mountain birch trees).

Soil sample collection took place on the 28^th and 29^th September 2019. Samples consisted of three to four soil cores (Ø 5 cm) from the upper 5 cm of the soil including both the litter and humus layers. After collection, the samples were separated into two halves that were stored at –18 °C. One half of the soil sample was used to determine soil water (SWC, drying at 105 °C for > 18 hours) and organic matter content (OM, loss on ignition at 405°C, 5 hr). The other half was freeze-dried, mixed manually, and a subset (2 ml) of this soil was homogenized before DNA extraction and C and N analyses. Soil C and N contents were analysed by vario Max CN Element Analyzer (Elementar; Langenselbold, Germany) and used to calculate soil C:N ratio and C and N stocks (SOC and SON, kg m^–2).

DNA extraction and quantitative PCR

Total DNA was extracted from 100 mg of dried soil using the NucleoSpin soil kit. Copy numbers of the fungal ITS2 region were quantified using the fungal primers gITS7, ITS4 and ITS4arch on the CFX Connect Real-Time System (Bio-Rad, CA, USA). Total ITS2 copy numbers were corrected for the proportion of non-fungal reads based on sequencing of the same region. Fungal abundance was then calculated as ITS2 copies g^–1 DW soil and ITS2 copies m^–2.

Soil fungal community and bioinformatics

The fungal ITS2 region was PCR-amplified using the same primers as above. Samples were sequenced at SciLifeLab NGI (Uppsala, Sweden) with PacBio Sequel I (Pacific Biosciences, Menlo park, CA, USA). The raw sequence data was quality filtered and clustered using the SCATA pipeline (https://scata.mykopat.slu.se). The final data set composed of 82,845 high-quality, non-unique fungal reads which clustered into 364 species-level operational taxonomical units (OTUs), and an average of 1,167 (222–8106) reads per sample. The OTUs were subjected to taxonomic identification using reference sequences from the UNITE database and from earlier annotated data sets (Clemmensen et al., 2015; Castaño, C., Lindahl, B., Clemmensen, K.E., unpublished data). The taxonomic affiliation and match to references associated with particular substrates was used to separate the OTUs into the following functional guilds: root-associated ascomycetes, root-associated basidiomycetes, moulds, yeasts, other saprotrophs (including litter-associated fungi), and lichenized fungi. The root-associated fungi were further classified into ectomycorrhizal fungi (EcM; including all fungal OTUs forming ectomycorrhizal symbiosis) and ericoid mycorrhizal fungi (ErM; excluding the OTUs that also form EcM). EcM fungi were further subdivided according to exploration type to cord-formers and simple mycelia as defined by Clemmensen et al. (2015). In total, 61% of the OTUs could be assigned to a functional guild, representing 86% of all reads.

Assessing fungal diversity

We used the package vegan (v2.5-5; Oksanen et al., 2019) in R version 3.6.3 (R Core Team, 2020) to assess alpha-diversity across samples. We used the read data rarefied to minimum sequence read number (n = 222) to calculate the following diversity indices: Species richness, Shannon–Wiener index (H), Simpson concentration (D), and Pielous’ evenness index. Shannon–Wiener and Simpson indices were further converted to indicate the effective number of species (i.e. Shannon-Wiener index is expressed as exponential function, exp(H), and Simpson index as the inverse, 1/D; Jost, 2006), and the effective Shannon–Wiener index was used in the calculation of Pielous’ evenness index.

Estimated plant abundance

Information on vegetation community composition in the study plots was derived from the plant reads captured by the ITS2 primers. 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 (Empetrum hermaphroditum, Phyllodoce caerulea and Calluna vulgaris) and bryophytes (genera Polytrichum and Dicranum). The four groups contributed to 89.0 ± 6.1% of plant cover in the blocks (mean ± SD; Ylänne, unpublished data). Metrics used in the calculation are found in the Data sheet 6: Read data on plant abundance).

Statistical analysis in R software

The code provided shows how 1) Diversity indices, 2) Grazing and tree effects on fungal diversity, plant and fungal abundance, and soil properties (linear mixed effects model, lme), 3) Grazing and tree effects on fungal community composition (multivariate generalized linear model, manyGLM), 4) Grazing and tree effects on fungal community composition (as deterred from non-metric multidimensional scaling of the fungal community, NMDS), 5) Correlations between soil fungal community and plant reads (linear models), and 6) Model comparison between soil fungal community and soil organic carbon stocks (comparison of linear models using Akaikes’ information criterion) were assessed.

All codes were run in R version 3.6.3 (R Core Team, 2020) using the packages vegan (v2.5-5; Oksanen et al., 2019), nlme (v3.1.144; Pinheiro et al., 2014), lsmeans (v2.30.0; Lenth, 2016), mvabund (v4.1.3; Wang et al., 2012), phyloseq (v1.30.0; McMurdie & Holmes, 2013), AICcmodavg (v2.2.2; Mazerolle, 2017) and, for data visualization, the package ggplot2 (v 3.3.0; Wickham, 2016).

References

Castaño C, Berlin A, Brandström Durling M, Ihrmark K, Lindahl BD, Stenlid J, Clemmensen KE, Olson Å (2020) Optimized metabarcoding with Pacific biosciences enables semi-quantitative analysis of fungal communities. New Phytologist 228: 1149–1158.

Clemmensen KE, Finlay RD, Dahlberg A, Stenlid J, Wardle DA, Lindahl BD (2015) Carbon sequestration is related to mycorrhizal fungal community shifts during long-term succession in boreal forests. New Phytologist 205: 1525–1536.

Jost L (2006) Entropy and diversity. Oikos 113: 363–375.

Lenth RV (2016) Least-Squares Means: The R Package lsmeans. Journal of Statistical Software 69: 1–33.

Mazerolle MJ (2019) AICcmodavg: Model Selection and Multimodel Inference Based on (Q)AIC(c). R package version 2.2-2.

McMurdie PJ, Holmes S (2013) Phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data. PLoS ONE 8: e61217.

Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P, et al. (2019) vegan: Community Ecology Package. R package version 2.5-2. Cran R.

Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team (2014) nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-118.

R Core Team (2020) R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing.

Wang Y, Naumann U, Wright ST, Warton DI (2012) Mvabund - an R package for model-based analysis of multivariate abundance data. Methods in Ecology and Evolution 3: 461–474.

Wickham H (2016) ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York.