https://orcid.org/0000-0002-8002-5122
Abstract
Land use and agricultural soil management affect soil fungal communities that ultimately influence soil health. Subsoils harbor nutrient reservoir for plants and can play a significant role in plant growth and soil carbon sequestration. Typically, microbial analyses are restricted to topsoil (0–30 cm) leaving subsoil fungal communities underexplored. To address this knowledge gap, we analyzed fungal communities in the vertical profile of four boreal soil treatments: long-term (24 years) organic and conventional crop rotation, meadow, and forest. Internal transcribed spacer (ITS2) amplicon sequencing revealed soil-layer-specific land use or agricultural soil management effects on fungal communities down to the deepest measured soil layer (40–80 cm). Compared to other treatments, higher proportion of symbiotrophs, saprotrophs, and pathotrophs + plant pathogens were found in forest, meadow and crop rotations, respectively. The proportion of arbuscular mycorrhizal fungi was higher in deeper (>20 cm) soil than in topsoil. Forest soil below 20 cm was dominated by fungal functional groups with proposed interactions with plants or other soil biota, whether symbiotrophic or pathotrophic. Ferrous oxide was an important factor shaping fungal communities throughout the vertical profile of meadow and cropping systems. Our results emphasize the importance of including subsoil in microbial community analyses in differently managed soils.
Introduction
Fungi play a key role in agricultural soil health by affecting soil structure through aggregation, nutrient cycling, plant health, and soil organic carbon (SOC) formation and decomposition (Powell and Rillig 2018, Toju et al. 2018, Bhattacharyya et al. 2022, Xiong and Lu 2022). High fungal diversity has been linked to improved soil multifunctionality and crop production (Wagg et al. 2011, Delgado-Baquerizo et al. 2016), and fungal abundance and activity to increased carbon sequestration into soil (Kallenbach et al. 2016, Bhattacharyya et al. 2022, Hannula and Morriën 2022). Soil management practices that promote diverse fungal communities in agricultural soil can potentially increase stable soil carbon formation and ultimately crop yield (Hannula and Morriën 2022), contributing to the UN sustainability goals for sustainable agriculture (United Nations 2015).
Fungi can be divided into functional groups according to their main mode of acquiring energy and nutrients (Nguyen et al. 2016). It has been proposed that, rather than an overall fungal community, certain fungal functional groups would better describe soil ecosystem functioning in agricultural soils (Ferris and Tuomisto 2015). Symbiotrophic fungi, especially arbuscular mycorrhizal fungi (AMF), which form interactions with plants and contribute to plant nutrient and water uptake (Smith and Read 2008), and saprotrophic fungi, which promote nutrient cycling in soil by decomposing organic material (Deacon et al. 2006), are potentially beneficial fungal groups for crop production. AMF have been linked to increased plant phosphorus uptake and ultimately higher plant productivity (van der Heijden et al. 1998, Fall et al. 2022) as well as to pathogen suppression in agricultural soils (Fall et al. 2022, Hannula and Morriën 2022). AMF can increase SOC by promoting plant photosynthate translocation into the soil matrix and by forming hyphal biomass (Jeewani et al. 2020, Parihar et al. 2020). Similarly, saprotrophs have been linked to higher soil fertility (Ning et al. 2021) and plant pathogen suppression (van der Wal et al. 2013). Saprotrophic fungi have been shown to increase SOC in forest ecosystems (Klink et al. 2022), and a similar effect in arable soils was recently proposed (Hannula and Morriën 2022). Investigating AMF and saprotrophs, as well as other fungal functional groups such as plant pathogens, that can have harmful effects on crop plant production (Corredor-Moreno and Saunders 2020), brings important knowledge on the health and functionality of agricultural soils. Agricultural management intensity, which is a measure of the fertilizer and biocide use, irrigation, and mechanization level (Foley et al. 2011), is known to affect soil microbial communities. For instance, high management intensity can decrease fungal biomass and the abundance of both AMF and saprotrophic fungi, likely due to their sensitivity to soil disturbance (Strickland and Rousk 2010, Thiele-Bruhn et al. 2012, Hydbom et al. 2017, Banerjee et al. 2019). Long-term organic management, representing a lowered management intensity in which chemical biocides and chemical fertilizers are not used, can promote fungal richness and abundance over more intensive conventional management (Martínez-García et al. 2018, Peltoniemi et al. 2021). Similarly, lowered management intensity in extensively managed grasslands or lands with permanent, predominantly herbaceous plant cover has been shown to promote fungal richness over arable lands (Banerjee et al. 2024). Yet, more knowledge of the effects of soil management on fungal communities across the soil vertical profile is needed as studies have mostly focused on topsoil, typically reaching to 20–30 cm depth at maximum (Henneron et al. 2022).
Agricultural soil carbon is decreasing globally (Lal et al. 2004). In Finland, agricultural mineral soils lose carbon at a yearly rate of 0.4% (Heikkinen et al. 2013). Globally, several agricultural management practices have been shown to promote SOC, including diverse crop rotation, organic amendments, no-tillage systems in some climate and soil type conditions, and organic farming, although the latter was recently revised to need additional actions, such as cover cropping and enhanced plant residue recycling (Francaviglia et al. 2017, Yang et al. 2019, Ogle et al. 2019, Zhang et al. 2021, Gaudaré et al. 2023). Soil management practices enabling the formation of extensive root systems and deep rooting plants can potentially increase SOC not only in topsoil but also in subsoil layers (below 30 cm), where half of the soil carbon of agricultural fields is stored (Balesdent et al. 2018, Hirte et al. 2021, Nguyen 2009, Paustian et al. 2016). In addition to roots, fungal hyphae contribute to the translocation of SOC deeper in the soil (Witzgall et al. 2021). The important role of deep soil layers in SOC sequestration (Button et al. 2022) and fungi in SOC dynamics further emphasizes the need to investigate fungal communities in the soil vertical profile to better understand the fate of SOC in agricultural soils.
Here we used amplicon sequencing of ribosomal RNA gene internal transcribed spacer (ITS2) to study fungal communities in the vertical profile of four soil treatments, organic and conventional cropping systems, unmanaged meadow, and forest, down to 80 cm deep after 24 years of field experiment. The comparison of organic and conventional treatments enables the assessment of the long-term combined effects of fertilizer type and herbicide usage on fungal communities. Organic treatment represents a less intensively managed system compared to conventional treatment, whereas meadow treatment represents the least intensive, natural-grassland-like management, creating a management intensity gradient from least intense to most intense: meadow–organic–conventional. Forest is included as a reference and represents the land use type that prevailed in the experiment area before conversion to agricultural system (Salonen et al. 2023), providing insight on how the transition into agricultural or meadow land use changes the fungal community over time. Our overall aim was to study how depth within the land use types (forest, meadow, and organic and conventional cropping systems) and soil management intensity within meadow–organic–conventional soil management intensity gradient influence fungal communities. In addition, we aimed to address which soil properties are the drivers of fungal community differences within the soil vertical profile. We hypothesize that lower soil management intensity in meadow and organic soils increase fungal diversity and promote the potentially beneficial AMF and saprotroph communities compared to more intensively managed conventional soil.
Materials and methods
Experimental site
The experimental site has been described in detail in Salonen et al. (2023). In short, the experimental site, established in 1995, located in south-western Finland (60°51′44.5′′N 23°31′24.5′′E) consists of (i) conventional cropping system (conventional); (ii) organic cropping system (organic); (iii) permanent meadow (meadow) all with three replicate field plots (Fig. S1). Additionally, an adjacent (iv) forest site (forest), which was sampled from one field plot only due to limited access of the soil coring machine, was included making up four different treatments. Soil in the plots was clayey and acidic, with a pH (measured in water) range of 5.29–6.06 in the topsoil, the forest plot being the most acidic. The organic and conventional plots had a 5-year crop rotation scheme with (i) a spring barley and a perennial grass mixture (timothy, meadow fescue, and red clover); (ii) a perennial grass mixture; (iii) a perennial grass mixture with winter rye sown in autumn; (iv) winter rye; and (v) a mixture of oat and garden pea (Salonen et al. 2023). The conventional plots were treated with herbicides approximately every other year and received mineral fertilizer, whereas the organic plots were fertilized with cattle slurry. Both the conventional and organic plots were periodically ploughed to 20 cm deep. The meadow plots were otherwise unmanaged, no fertilization, ploughing or cutting, only tree saplings were removed. The forest plot was primarily pine- and shrub-covered and had been regeneration felled 2 years prior to sampling.
Soil sampling
Soil samples were collected as explained in Salonen et al. (2023) in Autumn 2019. In brief, a motorized sample corer was used to insert cylinder tubes (100 cm long and 6.8 cm inner diameter) into the soil (Persson and Bergström 1991, Uusitalo et al. 2012). The vertical soil cores obtained were 60–80 cm long due to the inability of the tube to hold the deeper soil layers. Initially, three cores were sampled from each replicate plot, yielding nine soil cores (nine replicates) for each meadow, organic and conventional treatment, but for DNA analysis there was a reduced number of cores for meadow (eight cores) and organic (eight cores) treatments due to issues with sample storage. From forest, only three cores were taken from one plot. Immediately after sampling, cores were frozen at −20°C and later melted and cut into 10 cm-long fractions representing vertical soil layers and sieved through a 4 mm sieve. Mosses and leaves were removed from the top of forest cores before fractioning the soil layers. Soil for DNA analyses was kept at +4°C overnight before DNA extraction. Subsamples for various other soil analyses were taken and stored at 4°C or −20°C as required (Salonen et al. 2023).
Analysis of root biomass and soil chemical properties
Root biomass was measured from 300 g of wet soil by first separating roots by washing the soil in a 0.4-mm sieve and by manually removing nonroot organic material, after which roots were dried for 2 days at 60°C to constant weight and dry root weight per dry soil was calculated (g/kg). The analyses of soil chemical properties included: (i) pH (pH-H2O) measured in 1:2.5 (w:v) soil-deionized water-solution, (ii) nitrogen (N), and (iii) total carbon (C), which, under the relatively acidic soils in this study, represent organic C at least in the surface soil (Nelson and Sommers 1996), measured by dry combustion (LecoCHN 628, St. Joseph, Michigan, USA), (iv) C and N ratio (C/N) analyzed as C and N, (iv) dissolved organic carbon (DOC) measured from 1:10 (w:v) ratio of dried soil in deionized water filtered (0.2 μm, Nuclepore® polycarbonate, Whatman International Ltd) solution with a Shimadzu TOC-V CPH/CPN analyzer (Kyoto, Japan), (vi), iron and (vii) aluminum oxides (Fe-ox and Al-ox) extracted with acid ammonium oxalate extraction (Niskanen 1989), (viii) total phosphorus (P-tot) measured from aqua regia extraction (Crosland et al. 1995) facilitating microwave digestion with ICP-OES (inductively coupled plasma optical emission spectrometer, Thermo Scientific iCAP 6000 ICP spectrometer, Englanti), (ix) organic phosphorus (P-org) measured with a modified Saunders and Williams (1955) method (Kuo 1996) with a spectrophotometer (UVmini-1240 UV-VIS Spectrophotometer), (x) inorganic phosphorus (P-inorg) measured as the difference of P-tot and P-org, and (xi) water-extractable P (P-H2O) measured according to Hartikainen (1982) using 1:10 extraction ratio with discrete analyzer (Gallery Plus Discrete Analyzer, ThermoFisher Scientific, USA). All methods are described in more detail elsewhere (Salonen et al. 2023, Uhlgren 2023).
DNA extraction and sequencing
DNA was purified from 100 mg of soil using a NucleoSpin™ Soil Kit (Macherey & Nagel, Düren, Germany). Due to low DNA concentrations (Fig. S2), DNA from deep soil layers, 40–50 cm, 50–60 cm, 60–70 cm, and 70–80 cm, were pooled forming a fifth soil layer from 40 to 80 cm and resulting in a total of 140 samples (Table S1). The ITS2 DNA library of the samples was prepared in the Institute of Biotechnology, University of Helsinki. In short, DNA samples were first amplified with fITS7 and ITS4 primers (White et al. 1990, Ihrmark et al. 2012) that target the approximately 260 bp ITS2 region and then sequenced using paired end read chemistry on an Illumina MiSeq (Illumina, Inc., San Diego, CA, USA). Sequences were demultiplexed based on sample and the raw sequences were submitted to GenBank SRA (BioProject ID: PRJNA1131619).
Processing of sequence data
The Miseq reads were trimmed to remove adapters, and reads with unidentified bases or length shorter than 50 bp were discarded with Cutadapt software v3.4 (Martin 2011). Paired reads were joined with PEAR v0.9.8 (Zhang et al. 2014) and quality filtered with the VSEARCH v2.17.1 package (Rognes et al. 2016) to obtain 230–360 bp long reads with a maximum expected error of 1 counted from phred quality score. Sequences were dereplicated and de novo chimera detection was performed using the USEARCH v11 UCHIME algorithm (Edgar et al. 2011). The sequences were clustered into operational taxonomic units (OTUs) with VSEARCH using ≥0.97 pairwise sequence identity with the cluster centroid. The quality-filtered reads were mapped back to the OTUs with ≥0.97 identity to obtain a count table. OTUs were taxonomically classified against the UNITEv10_sh_99 database (Nilsson et al. 2019) using Mothur v1.45.3 software (Schloss et al. 2009) and 60% bootstrap confidence to obtain a taxonomy table. Singletons and OTUs classified as nonfungal or as unclassified at phylum level were filtered out from the count and taxonomy tables. OTU counts were normalized by library size (relative abundance transformation).
Two additional classifications of the OTUs were done. First, to study AMF communities, OTUs belonging to the Glomeromycota phylum were assigned to an AMF guild, edaphophilic, rhizophilic, or ancestral, based on family level classification in Weber et al. (2019) and classifying Entrophosporaceae according to the former family Claroideoglomeraceae (Błaszkowski et al. 2022). Second, to analyze fungal functional groups, fungal taxa were assigned to functional trophic modes and functional guilds with the FUNGuild database (Nguyen et al. 2016; downloaded 30 May 2024), including taxa with a confidence ranking level "highly probable" or "probable." The majority of fungal OTUs (55.2%) and reads (69.5%) were assigned to a functional group so that more reads were assigned in meadow (76.2%) and forest (71.7%) than in organic (68.2%) and conventional (64.1%) samples. Unimodal (i.e. having only one specified guild) ectomycorrhizal and endophyte guilds, and unimodal and multimodal arbuscular mycorrhizal guilds (i.e. arbuscular mycorrhizal and arbuscular mycorrhizal–endophyte) and unimodal plant pathogen guild were specified and are referred from here on as ectomycorrhizal, endophyte, arbuscular mycorrhizal and plant pathogen groups.
Statistical analysis
All analyses of the fungal OTUs were done in RStudio with R version 4.2.2 (R Core Team 2022). RStudio software code on which the conclusions of the paper rely is available on GitHub (https://github.com/lehakkin/Fungal-ITS-soil-vertical-profile). Dissimilarities in fungal communities were analyzed by calculating the Bray–Curtis distance of samples with vegdist-function from vegan package v2.6–2 (Oksanen et al. 2012) using OTU relative abundances. The obtained distance matrix was used for principal coordinates analysis (PCoA) with function cmdscale (R Core Team 2022) using eigen-values and three dimensions (k = 3). To visualize beta diversity, an ordination plot was produced with vegan's ordiplot-function (Oksanen et al. 2012). Vectors of soil properties, C, N, C/N, DOC, pH, total P-tot, P-inorg, P-org, P-H2O, Fe-ox, Al-ox, and root biomass, were fitted to the ordination data with Vegan's envfit-function with 999 permutations (Oksanen et al. 2012) and only significant (P < .05) properties were included in the figure. Species scores of the 0.1% most abundant and 100% of the best-fitting fungal OTUs were picked with goeveg's ordiselect-function (von Lampe and Schellenberg 2023).
To partition out the effect of soil layer and treatment on fungal communities, a permutational multivariate analysis of variance (PERMANOVA) was carried out with Bray–Curtis distance using vegan's adonis2-function (Oksanen et al. 2012) with 9999 permutations. To investigate pairwise differences between treatments and soil layers, pairwise PERAMANOVA was performed. For this, treatment and soil layer were checked for homogeneity of variance with betadisper-function followed by one-way analysis of variance (ANOVA) and Tukey tests. Variances were similar between all treatments (ANOVA; P = .409) and soil layer (Tukey; P > .1) except for the deepest soil layer, 40–80 cm, which had a significantly different variance compared to all other soil layers (Tukey; P < .05). Additionally, PERMANOVA was used to investigate the effect of each soil property on fungal communities in the soil vertical profile in meadow, organic, and conventional treatments. Compared to some other tests (e.g. Mantel's test), PERMANOVA deals better with the heterogeneity of dispersion (Anderson and Walsh 2013) that was present in the deepest soil layer.
To identify differentially abundant fungal genera in meadow, organic, and conventional treatments, first the 20 most abundant genera in each of the five soil layers of the three treatments (15 groups) were chosen, adding up to 96 unique genera. Second, differently abundant fungal genera between relevant comparisons (successive soil layers within treatment and the same soil layers between treatments) were identified with Wilcoxon test of the relative abundances (Table S3), adding up to 70 differently abundant fungal genera that were used to construct a heatmap with z-scores of the relative abundances using package pheatmap v1.0.12 (Kolde 2019). For the forest treatment, all the 10 most abundant fungal genera in each forest soil layer, adding up to 31 genera, were used to construct a heatmap with z-scores of the relative abundances.
Fungal and fungal functional group richness (observed richness, i.e. number of OTUs) were calculated with package microbiome v1.18.0 (Lahti et al. 2017). Correlation of richness and fungal taxa proportion with soil properties and depth was assessed with Spearman's rank test (Spearman). Differences in fungal taxa proportion and fungal richness between treatments and soil layers were conducted by first examining homogeneity of variances with Levene's test and in case of homoscedasticity, ANOVA test was conducted followed by a Tukey Honest Significant Difference test (Tukey) for multiple pairwise comparison between the means of groups. In case of heteroscedasticity, a nonparametric Kruskal–Wallis ANOVA test (Kruskal–Wallis) was conducted, followed by a pairwise Wilcoxon test (Wilcoxon). The statistical significance level was set to P ≤ .05 in all analyses, and for adjusted P values the Benjamini–Hochberg method was used except in pairwise PERMANOVA where Bonferroni adjustment was used instead.
Results
Fungal community differences across soil profile and land use or agricultural management types
We obtained 11 539 503 filtered fungal reads without singletons which clustered into 20 610 OTUs in the 140 samples. Apart from the deepest soil layer (40–80 cm; 36 000 reads mean), all soil layers had a similar mean number of reads per sample (81 000–101 000 reads mean; Table S2). The smallest number of reads, 44, was found in one conventional treatment 40–80 cm sample, but it was kept in the analysis due to having similar community composition and richness values as its counterpart samples (data not shown). Otherwise, reads per sample were 526 or more. A small proportion of OTUs was found across all treatments (12.5%), cropping systems (organic and conventional; 22.1%), or soil layers (4.9%) (Fig. 1). OTUs unique to each soil layer decreased toward deeper soil layers so that majority of the OTUs in the dataset (97.7%) were found in the first three soil layers. Fungal community composition was affected by the soil layer (R2 = 0.18; P < .001), treatment (R2 = 0.10; P < .001), and their interaction (R2 = 0.13; P < .001) as indicated by PCoA (Fig. 2; Fig. S3) and confirmed by PERMANOVA (Table S4A). Each soil layer in the meadow differed from comparable layers in conventional and organic treatments in soil layers between 0 and 40 cm. Organic and conventional treatments differed in the first two and the deepest soil layers (0–10 cm; 10–20 cm; 40–80 cm) (Table S4B–F). Forest was not included in the pairwise soil layer-specific analysis due to lower replicate numbers (3). Within each treatment, fungal communities differed between most of the soil layers, but less differences were found between successive soil layers (Table S4G–I). Successive soil layers differed the most in meadow followed by organic treatments, and no differences were found in conventional treatment.
Venn diagram of shared and unique fungal OTUs in forest, meadow, organic, and conventional treatments (left) and in the five different soil layers (right).
Principal coordinates analysis (PCoA) ordination with Bray–Curtis distance showing the first and second axes (PC1 and PC2). Each data point is a sample. Colors represent the four treatments and the shapes of the five soil layers. Significant environmental variables fitted with vegan's envfit-package and the 0.1% most abundant best fitting fungal OTUs based on species scores picked with ordiselect-function are shown in the figure. PCoA with first and third axes is presented in the Fig. S3.
We used PERMANOVA to assess the effect of individual soil properties on fungal communities within each soil layer and in the whole soil profile irrespective of treatment (Table 1). Only the more similar meadow, organic, and conventional were included in the analysis. Soil properties affected fungal communities in a soil-layer-specific manner, and most strongly in the first two soil layers (0–10 cm and 10–20 cm). In the whole soil profile (0–80 cm), the strongest effects (R2 ≥ 0.13) were by P-org, N, C, and pH, whereas P-tot, root biomass, DOC, C/N, and Fe-ox had an intermediate effect (R2 = 0.10–0.12), and P-H2O and P-inorg had small effect (R2 = 0.01–0.05). However, the only soil properties consistently influencing fungal communities between 0 and 40 cm when soil layers were analyzed individually were root biomass, pH, and Fe-ox, out of which Fe-ox continued to have an effect also in the last soil layer (40–80 cm).
PERMANOVA R2 values and significance levels with 9999 permutations between soil properties and fungal Bray–Curtis distance in each soil layer separately and in the whole soil profile of meadow, organic and conventional treatment.
| Soil property | 0–10 cm | 10–20 cm | 20–30 cm | 30–40 cm | 40–80 cm | 0–80 cm |
|---|---|---|---|---|---|---|
| P-org | 0.158** | 0.054(ns) | 0.090** | 0.079* | 0.039(ns) | 0.145*** |
| N | 0.272*** | 0.192** | 0.097** | 0.044(ns) | 0.038(ns) | 0.144*** |
| C | 0.281*** | 0.240*** | 0.100** | 0.043(ns) | 0.047(ns) | 0.140*** |
| pH | 0.110** | 0.182*** | 0.118*** | 0.071* | 0.051(ns) | 0.131*** |
| P-tot | 0.046(ns) | 0.049(ns) | 0.065(ns) | 0.055(ns) | 0.031(ns) | 0.120*** |
| Log10 root biomass | 0.159*** | 0.141** | 0.103** | 0.078** | 0.040(ns) | 0.116*** |
| DOC | 0.258*** | 0.215*** | 0.092** | 0.053(ns) | 0.030(ns) | 0.115*** |
| C/N | 0.218*** | 0.113** | 0.063(ns) | 0.063(ns) | 0.051(ns) | 0.106*** |
| Fe-ox | 0.139** | 0.225*** | 0.107** | 0.079* | 0.080** | 0.105*** |
| P-H2O | 0.092* | 0.070(ns) | 0.069* | 0.041(ns) | 0.046(ns) | 0.042*** |
| P-inorg | 0.068(ns) | 0.076* | 0.042(ns) | 0.039(ns) | 0.029(ns) | 0.015* |
| All | 0.338*** |
| Soil property | 0–10 cm | 10–20 cm | 20–30 cm | 30–40 cm | 40–80 cm | 0–80 cm |
|---|---|---|---|---|---|---|
| P-org | 0.158** | 0.054(ns) | 0.090** | 0.079* | 0.039(ns) | 0.145*** |
| N | 0.272*** | 0.192** | 0.097** | 0.044(ns) | 0.038(ns) | 0.144*** |
| C | 0.281*** | 0.240*** | 0.100** | 0.043(ns) | 0.047(ns) | 0.140*** |
| pH | 0.110** | 0.182*** | 0.118*** | 0.071* | 0.051(ns) | 0.131*** |
| P-tot | 0.046(ns) | 0.049(ns) | 0.065(ns) | 0.055(ns) | 0.031(ns) | 0.120*** |
| Log10 root biomass | 0.159*** | 0.141** | 0.103** | 0.078** | 0.040(ns) | 0.116*** |
| DOC | 0.258*** | 0.215*** | 0.092** | 0.053(ns) | 0.030(ns) | 0.115*** |
| C/N | 0.218*** | 0.113** | 0.063(ns) | 0.063(ns) | 0.051(ns) | 0.106*** |
| Fe-ox | 0.139** | 0.225*** | 0.107** | 0.079* | 0.080** | 0.105*** |
| P-H2O | 0.092* | 0.070(ns) | 0.069* | 0.041(ns) | 0.046(ns) | 0.042*** |
| P-inorg | 0.068(ns) | 0.076* | 0.042(ns) | 0.039(ns) | 0.029(ns) | 0.015* |
| All | 0.338*** |
For the whole soil profile, the R2 value of a PERMANOVA model including all soil properties is shown at the bottom line (“All”). Significance codes: *** < .001 ** < .01 * < .05.
PERMANOVA R2 values and significance levels with 9999 permutations between soil properties and fungal Bray–Curtis distance in each soil layer separately and in the whole soil profile of meadow, organic and conventional treatment.
| Soil property | 0–10 cm | 10–20 cm | 20–30 cm | 30–40 cm | 40–80 cm | 0–80 cm |
|---|---|---|---|---|---|---|
| P-org | 0.158** | 0.054(ns) | 0.090** | 0.079* | 0.039(ns) | 0.145*** |
| N | 0.272*** | 0.192** | 0.097** | 0.044(ns) | 0.038(ns) | 0.144*** |
| C | 0.281*** | 0.240*** | 0.100** | 0.043(ns) | 0.047(ns) | 0.140*** |
| pH | 0.110** | 0.182*** | 0.118*** | 0.071* | 0.051(ns) | 0.131*** |
| P-tot | 0.046(ns) | 0.049(ns) | 0.065(ns) | 0.055(ns) | 0.031(ns) | 0.120*** |
| Log10 root biomass | 0.159*** | 0.141** | 0.103** | 0.078** | 0.040(ns) | 0.116*** |
| DOC | 0.258*** | 0.215*** | 0.092** | 0.053(ns) | 0.030(ns) | 0.115*** |
| C/N | 0.218*** | 0.113** | 0.063(ns) | 0.063(ns) | 0.051(ns) | 0.106*** |
| Fe-ox | 0.139** | 0.225*** | 0.107** | 0.079* | 0.080** | 0.105*** |
| P-H2O | 0.092* | 0.070(ns) | 0.069* | 0.041(ns) | 0.046(ns) | 0.042*** |
| P-inorg | 0.068(ns) | 0.076* | 0.042(ns) | 0.039(ns) | 0.029(ns) | 0.015* |
| All | 0.338*** |
| Soil property | 0–10 cm | 10–20 cm | 20–30 cm | 30–40 cm | 40–80 cm | 0–80 cm |
|---|---|---|---|---|---|---|
| P-org | 0.158** | 0.054(ns) | 0.090** | 0.079* | 0.039(ns) | 0.145*** |
| N | 0.272*** | 0.192** | 0.097** | 0.044(ns) | 0.038(ns) | 0.144*** |
| C | 0.281*** | 0.240*** | 0.100** | 0.043(ns) | 0.047(ns) | 0.140*** |
| pH | 0.110** | 0.182*** | 0.118*** | 0.071* | 0.051(ns) | 0.131*** |
| P-tot | 0.046(ns) | 0.049(ns) | 0.065(ns) | 0.055(ns) | 0.031(ns) | 0.120*** |
| Log10 root biomass | 0.159*** | 0.141** | 0.103** | 0.078** | 0.040(ns) | 0.116*** |
| DOC | 0.258*** | 0.215*** | 0.092** | 0.053(ns) | 0.030(ns) | 0.115*** |
| C/N | 0.218*** | 0.113** | 0.063(ns) | 0.063(ns) | 0.051(ns) | 0.106*** |
| Fe-ox | 0.139** | 0.225*** | 0.107** | 0.079* | 0.080** | 0.105*** |
| P-H2O | 0.092* | 0.070(ns) | 0.069* | 0.041(ns) | 0.046(ns) | 0.042*** |
| P-inorg | 0.068(ns) | 0.076* | 0.042(ns) | 0.039(ns) | 0.029(ns) | 0.015* |
| All | 0.338*** |
For the whole soil profile, the R2 value of a PERMANOVA model including all soil properties is shown at the bottom line (“All”). Significance codes: *** < .001 ** < .01 * < .05.
Taxonomic composition of fungi
The proportion of the five most abundant fungal phyla and classes differed the most between forest and other treatments, but differences were also found between meadow, organic, and conventional treatments (Fig. 3A and B). Within the whole soil profile, forest was dominated by Basidiomycota (53%), and meadow, organic, and conventional by Ascomycota (67%–71%) whereas the proportion of Mortierellomycota, Glomeromycota, and Rozellomycota did not differ significantly between treatments (Wilcoxon P > .05; Fig. 3A). Soil layer did not significantly affect the proportions of phyla in the forest (Wilcoxon; P > .05; Fig. 3A; Table S10). However, in meadow, organic, and conventional treatments, the proportion of Basidiomycota was higher in the first three soil layers (0–30 cm) compared to the two deepest soil layers (30–80 cm), Mortierellomycota proportion peaked at the fourth soil layer (30–40 cm) being higher than in any other soil layer, Glomeromycota proportion grew toward deeper soil layers within the 0–40 cm soil profile, Chytridiomycota proportion was low throughout the soil profile but higher in first three soil layers compared to the two deepest soil layers (Wilcoxon; P < .05; Table S9), and Ascomycota did not differ between soil layers (Wilcoxon; P > .05).
Proportion of fungal taxa and functional groups. Composition bar plots indicating the proportions of (A) fungal phyla, (B) class, and (C) functional groups in the five soil layers of the four treatments. In (A) and (B), fungal phyla or classes with low relative abundance are grouped in “Other.” In (C) ectomycorrhizal, arbuscular mycorrhizal, and endophyte groups within the symbiotroph trophic mode, and plant pathogen group within the pathotroph trophic mode are specified. (D) Heatmap of most abundant significantly differently abundant 70 fungal genera in meadow, organic, and conventional treatment soil layers. (E) Heatmap of most abundant fungal genera in forest soil layers. The color of each cell ranging from blue to white to red indicates the z-score of the mean proportion in each treatment layer. Forest and other treatments were kept separately in heatmaps to better visualize differences between the more similar treatments: meadow, organic, and conventional.
Among common fungal classes, the proportion of Leotiomycetes and Tremellomycetes did not differ between treatments (Kruskal; P > .117), whereas Sordariomycetes was higher in meadow, organic, and conventional compared to forest and higher in organic and conventional compared to meadow, Dothideomycetes was higher in meadow and conventional compared to forest, and Agaricomycetes was higher in forest than in other treatments and also higher in organic compared to conventional treatment (Wilcoxon; P < 0.05; Fig. 3B). In the soil vertical profile of meadow, organic, and conventional treatments, Leotiomycetes was more associated with below 20 cm (20–80 cm) than the first two soil layers (0–20 cm), Sordariomycetes with the first three soil layers (0–30 cm) than fourth soil layer (30–40 cm), Dothideomycetes and Tremellomycetes with first three (0–30 cm) than deeper soil layers (20–80 cm), and Mortierellomycetes with third and fourth (20–40 cm) than other soil layers (Wilcoxon; P < .05; Table S11). Soil layer did not significantly affect the proportions of most abundant fungal classes in the forest (Wilcoxon; P > .05; Table S12). Fungal classes that mostly had low proportion but had higher proportion in specific layers of specific treatment were Pezizomycetes in forest 30–40 cm (ANOVA; P = .006), Geoglossomycetes in meadow 10–40 cm (Kruskal; P = .000), Microbotryomycetes in conventional 40–80 cm (Kruskal, P = .030), and Orbiliomycetes in meadow 10–30 cm (Kruskal; P = .000).
The proportions of most abundant fungal genera that differed significantly in relevant treatment or soil layer comparisons were visualized in the heatmap (Fig. 3D). Heatmap indicated distinctive genera for organic treatment including Cladorrhinum, Cheilymenia, Paraphoma, Apiotrichum, Epicoccum, Enterocarpus, and Podospora, distinctive genera to conventional treatment type including Exophiala, Fusicolla, Trichoderma, Linnemannia, Gibellulopsis, and Apodus, and distinctive genera to meadow including Leohumicola, Paraphaeosphaeria, Serendipita, Mycena, Laburnicola, and Tolypocladium. In forest, the proportions of most of the abundant fungal genera showed soil layer association trends (Fig. 3E) but these could not be validated with statistical tests (data not shown).
Correlations between depth and proportion of the five most abundant phyla, classes, and genera are shown in Table S5.
Fungal functional groups in treatments across soil vertical profile
At the fungal trophic mode level, forest was associated with higher symbiotroph than other treatments and low saprotroph proportion (Fig. 3C; Table S6). In meadow, organic, and conventional treatments, saprotrophs were the dominating group throughout the soil profile, and they were more associated with meadow than organic and conventional treatments. Organic and conventional treatments had higher pathotroph, plant pathogen, and endophyte proportions than forest and meadow. In the soil vertical profile of meadow, organic, and conventional treatments, symbiotroph proportion was low until the two deepest soil layers (30–40 cm and 40–80 cm) and was positively affected by depth (Spearman; rho = 0.257; P = .000), whereas pathotroph proportion was higher in topsoil layers and was strongly negatively affected by depth (Spearman; rho = −0.660; P = .000; Table S5). In addition, the saprotroph proportion was negatively affected by depth (Spearman; rho = −0.455; P = .000) but the lowest proportion was in the 30–40 cm soil layer (16.4%) rather than in the deepest soil layer (40–80 cm) (24.5%).
The main symbiotrophic group in organic, conventional, and meadow soil, AMF belonging to the Glomeromycota phylum was investigated in more detail. At the genus level, AMF communities differed between soil layers (PERMANOVA; R2 = 0.084; P = 0.000) and between treatments (PERMANOVA; R2 = 0.114; P = 0.000). However, in pairwise comparison, forest (R2 = 0.143; R2 = 0.070) and meadow (R2 = 0.136; R2 = 0.056) differed from organic and conventional, but organic could not be separated from conventional or forest from meadow (Table S13). Significant differences in the proportion of AMF genus, family, order, or guild between treatments were not found (Wilcoxon; P > .05) other than Ambisporaceae, which was higher in forest compared to other treatments (Wilcoxon; P < .05) (Fig. 4).
The proportions of different AMF taxa belonging to Glomeromycota phylum out of all fungi in the five soil layers of the four treatments. In the bubble chart, AMF genera or higher taxonomic level is indicated on the y axis, soil layer in the x axis and treatments are separated by facets. The color of the bubble indicates fungal order or higher taxonomic rank if order level identification was not available, and the size of the bubble indicates proportion. Font color in the y axis indicates AMF guild. Characters in parentheses: phylum (p), family (f), order (o), or genus (g).
Fungal alpha diversity
Fungal alpha diversity assessed as richness was highest at the two first soil layers (0–20 cm), decreasing toward deeper soil layers, except in meadow where fungal richness was highest in the second and third soil layers (10–30 cm) (Fig. 5A; Table S7). Between meadow, organic, and conventional treatment, fungal richness differed significantly only in the first soil layer (0–10 cm) where organic and conventional treatments had higher richness compared to meadow, and in the deepest soil layer (40–80 cm) where organic had higher richness compared to conventional (Tukey; P < .05). Between meadow, organic, and conventional, the AMF richness differed significantly in 20–30 cm and 30–40 cm soil layers where meadow was higher compared to conventional treatment (Tukey; P < .05), saprotroph richness differed in the 0–10 cm soil layer where organic was higher compared to meadow (Tukey; P < .05), and pathotroph richness differed in the first two soil layers (0–10 cm and 10–20 cm) where organic and conventional treatments were higher compared to meadow (Tukey; P < .05).
Diversity of fungi and selected fungal functional groups in the five soil layers of the four treatments. In (A), the overall fungal richness is given as thousands of OTUs. Richness of selected functional groups based on FUNGuild annotation (Nguyen et al. 2016), (B) richness of arbuscular mycorrhizal fungi (AMF), (C) saprotrophic fungi, and (D) pathotrophic fungi. Mean values and standard error of mean bars are shown. The significance of the statistical comparison between treatments within each soil layer (Tukey Honest Significant Difference test) is given with lower case letters so that shared letter denotes no difference.
Spearman's rank test with meadow, organic, and conventional samples indicated fungal richness correlating positively the strongest (rho > 0.5) with C, N, DOC, P-org, P-tot, C/N, root biomass, and Fe-ox, and negatively the strongest (rho < –0.5) with depth and pH (Table 2). AMF richness correlated with the same soil properties as fungal richness, but weaker. The only exception was Fe-ox, which correlated more strongly with AMF richness than with fungal richness (rho = 0.676; rho = 0.552). AMF relative abundance was affected differently by soil properties than AMF richness, so that apart from the positive effect of depth, only different forms of P, including P-org, P-tot, and P-H2O, but not P-inorg, had a significant and negative effect on AMF relative abundance.
Spearman's rank correlation rho values and significance levels between soil properties and fungal richness, AMF richness, and AMF relative abundance (AMF RA) in meadow, organic, and conventional treatments.
| Spearman's rank correlation rho | |||
|---|---|---|---|
| Soil property | Fungal richness | AMF richness | AMF RA |
| C | 0.696*** | 0.510*** | −0.135(ns) |
| N | 0.682*** | 0.458*** | −0.180* |
| DOC | 0.671*** | 0.516*** | −0.158(ns) |
| P-org | 0.659*** | 0.357*** | −0.242** |
| pH | −0.670*** | −0.574*** | 0.062(ns) |
| P-tot | 0.587*** | 0.224* | −0.324*** |
| C/N | 0.624*** | 0.562*** | −0.041(ns) |
| Log10 root biomass | 0.548*** | 0.417*** | −0.133(ns) |
| Fe-ox | 0.553*** | 0.619*** | 0.016(ns) |
| P-H2O | 0.188* | −0.010(ns) | −0.395*** |
| P-inorg | −0.079(ns) | −0.181* | −0.072(ns) |
| depth | −0.776*** | −0.393*** | 0.209* |
| Spearman's rank correlation rho | |||
|---|---|---|---|
| Soil property | Fungal richness | AMF richness | AMF RA |
| C | 0.696*** | 0.510*** | −0.135(ns) |
| N | 0.682*** | 0.458*** | −0.180* |
| DOC | 0.671*** | 0.516*** | −0.158(ns) |
| P-org | 0.659*** | 0.357*** | −0.242** |
| pH | −0.670*** | −0.574*** | 0.062(ns) |
| P-tot | 0.587*** | 0.224* | −0.324*** |
| C/N | 0.624*** | 0.562*** | −0.041(ns) |
| Log10 root biomass | 0.548*** | 0.417*** | −0.133(ns) |
| Fe-ox | 0.553*** | 0.619*** | 0.016(ns) |
| P-H2O | 0.188* | −0.010(ns) | −0.395*** |
| P-inorg | −0.079(ns) | −0.181* | −0.072(ns) |
| depth | −0.776*** | −0.393*** | 0.209* |
Significance codes: *** < .001 ** < .01 * < .05.
Spearman's rank correlation rho values and significance levels between soil properties and fungal richness, AMF richness, and AMF relative abundance (AMF RA) in meadow, organic, and conventional treatments.
| Spearman's rank correlation rho | |||
|---|---|---|---|
| Soil property | Fungal richness | AMF richness | AMF RA |
| C | 0.696*** | 0.510*** | −0.135(ns) |
| N | 0.682*** | 0.458*** | −0.180* |
| DOC | 0.671*** | 0.516*** | −0.158(ns) |
| P-org | 0.659*** | 0.357*** | −0.242** |
| pH | −0.670*** | −0.574*** | 0.062(ns) |
| P-tot | 0.587*** | 0.224* | −0.324*** |
| C/N | 0.624*** | 0.562*** | −0.041(ns) |
| Log10 root biomass | 0.548*** | 0.417*** | −0.133(ns) |
| Fe-ox | 0.553*** | 0.619*** | 0.016(ns) |
| P-H2O | 0.188* | −0.010(ns) | −0.395*** |
| P-inorg | −0.079(ns) | −0.181* | −0.072(ns) |
| depth | −0.776*** | −0.393*** | 0.209* |
| Spearman's rank correlation rho | |||
|---|---|---|---|
| Soil property | Fungal richness | AMF richness | AMF RA |
| C | 0.696*** | 0.510*** | −0.135(ns) |
| N | 0.682*** | 0.458*** | −0.180* |
| DOC | 0.671*** | 0.516*** | −0.158(ns) |
| P-org | 0.659*** | 0.357*** | −0.242** |
| pH | −0.670*** | −0.574*** | 0.062(ns) |
| P-tot | 0.587*** | 0.224* | −0.324*** |
| C/N | 0.624*** | 0.562*** | −0.041(ns) |
| Log10 root biomass | 0.548*** | 0.417*** | −0.133(ns) |
| Fe-ox | 0.553*** | 0.619*** | 0.016(ns) |
| P-H2O | 0.188* | −0.010(ns) | −0.395*** |
| P-inorg | −0.079(ns) | −0.181* | −0.072(ns) |
| depth | −0.776*** | −0.393*** | 0.209* |
Significance codes: *** < .001 ** < .01 * < .05.
Discussion
In a recent meta-analysis, it was shown that in the deeper soil layers, there are on average 47% of soil organic C stocks of agricultural fields (Balesdent et al. 2018). Similarly in the forest soils, it has been shown that the total soil C stock under 20 cm may be up to 50% of the total (Jobbágy and Jackson 2000) and up to 75% of SOM can be found in subsoil (B and C horizons) (Rumpel et al. 2002). Considering deeper soil layers as reservoirs for C, different agricultural management practices can have a significant role both as enhancing fresh C input into deeper layers (Lessmann et al. 2022, Gaudaré et al. 2023), as well as modifying the microbial communities responsible for SOC decomposition and plant nutrient uptake (Morugan-Coronado et al. 2022). However, we still lack a comprehensive view of how land use or soil management influences microbial communities in the soil vertical profile and how this ultimately affects the fate of SOC.
Depth together with land use and agricultural soil management affected fungal community composition
The analysis of the vertical soil profile of the four treatments showed that fungal communities were affected by soil layer and treatment and the treatment effect varied between the studied five soil layers. Overall, we found soil layer to have a bigger effect on fungal community differences compared to treatment. Fungal community composition and diversity have previously been shown to be influenced by depth in cropping systems (Schlatter et al. 2018, Yin et al. 2021) and forest (Baldrian et al. 2012). Similarly, there are numerous studies showing how agricultural soil management intensity shapes fungal communities in topsoil (Sun et al. 2016, Gottshall et al. 2017, Vahter et al. 2022, Wu et al. 2022). However, previously the comparison of organic and conventional treatment effects on the fungal community has been done down to 30 cm, but we lack studies where below 30 cm layers are analysed (Epp Schmidt et al. 2022). Here, we show that the treatment effect between organic and conventional cropping systems can be seen down to the deepest measured soil layer 40–80 cm (Table S4B–F). Conventional and organic plots had the same 5-year crop rotation and three different crops growing during the sampling year, indicating that agricultural management affects fungal communities regardless of the crop.
Fungal richness was not negatively associated with soil management intensity
In line with a previous study by Schlatter et al. (2018), our results on fungal richness showed a consistent decrease in relation to depth in soil layers between 10–80 cm in all treatments. Fungal richness in meadow, organic, and conventional treatments differed in topsoil (0–10 cm) where organic and conventional had more diverse fungal community compared to meadow and in the deepest soil layer (40–80 cm) where organic treatment had more diverse fungal community compared to conventional. Interestingly, contrary to what we hypothesized and what has been found in multiple previous studies (Martínez-García et al. 2018, Peltoniemi et al. 2021, Banerjee et al. 2024), low management intensity did not promote higher fungal richness in topsoil. This, however, follows the somewhat surprising fungal diversity pattern found in a Europe-wide study across land-use intensity gradients (woodland–grassland–cropland), where higher land use intensity correlated with higher fungal diversity (Labouyrie et al. 2023). Similarly, in grasslands, the intensification of land management practices has been found to have either neutral or positive effects on belowground fungal diversity (Allan et al. 2014, Gossner et al. 2016). In diverse environments such as the meadow, organic, and conventional soils in our study, the common understanding in ecology that a higher species richness contributes to higher ecosystem functioning (Loreau et al. 2001) has been disputed (Nielsen et al. 2011). Ecosystem functions have rather been linked to succession of fungal communities than to high OTU richness (Hoppe et al. 2016). We do not have data for temporal succession in our soils, but we know that fungal communities were more specialized vertically in meadow (Table S4G), probably due to higher litter input and the lack of interruption by periodic ploughing. This spatial specialization in meadow could possibly lower fungal diversity in individual soil layers. In addition, the lower topsoil pH in the 0–10 cm soil layer of meadow compared to organic treatment and marginally compared to conventional management may have attributed to the lower fungal diversity in meadow (Zheng et al. 2019). The over two-fold higher DOC in the 0–10 cm soil layer of meadow, most probably caused by the high litter input, may also have lowered topsoil fungal richness in meadow similarly as in a previous study where higher arable soil DOC and lower fungal richness were found in straw mulch soil compared to soil without mulch (Huang et al. 2019). We did not find difference in fungal richness between organic and conventional in the first four soil layers (0–40 cm). Similarly, in a study with organically fertilized (pig manure) and chemically fertilized crop field, and in long-term organic and conventional cereal crop systems, no significant differences in fungal Shannon diversity (Suleiman et al. 2019) or OTU richness (Peltoniemi et al. 2021) between the management types were found, but rather in the fungal ITS2 copy numbers (Peltoniemi et al. 2021), indicating that management effect on fungal diversity could be more subtle compared to fungal abundance, which was not measured in this study. However, our study provides only a single time point view of fungal diversity which can change during the growing season and between years (Degrune et al. 2017). Considering our findings and the literature, the overall effect of management intensity on fungal richness remains somewhat unclear.
Management intensity affected AMF richness below the surface soil
Previously, it has been shown that rather than the overall fungal community, specialized microbial groups are linked to soil ecosystem functioning and may better describe the effects of land use or soil management intensity (Wang et al. 2022). Symbiotrophic fungi in general and specifically AMF can benefit plant productivity and soil fertility (van der Heijden et al. 1998, Smith and Read 2008 et al. 2008, Jeewani et al. 2020, Parihar et al. 2020, Fall et al. 2022, Hannula and Morriën 2022). Although lower agricultural soil management intensity is shown to positively affect AMF (Hydbom et al. 2017, Banerjee et al. 2019), we did not find a significant effect of treatment on AMF or symbiotroph proportion between meadow, organic, and conventional treatments. AMF richness, however, differed between the low-intensity meadow and the highest intensity conventional treatment in the 20–30 and 30–40 cm soil layers. Organic soil which represents a lowered management intensity fell between the intensity extremes and could not be statistically differentiated from either.
The management intensity effect on AMF richness can be attributed to different management practices. For instance, AMF are shown to be negatively affected by fertilization overall (Hannula et al. 2021, Luo et al. 2021), and the use of mineral fertilization over manure can further suppress AMF (Wang et al. 2018), which could explain higher AMF richness in unfertilized meadow compared to mineral-fertilized conventional treatment. The differences in root biomass between treatments which followed the management intensity gradient (higher root biomass in lower management intensity; Fig. S5; Table S8) and the lack of disturbance related to tillage operations in meadow may have promoted higher AMF richness in meadow (Hiiesalu et al. 2014, Schmidt et al. 2019). Plant diversity was not measured from the treatment plots in the sampling year (2019), so we cannot fully assess the effect of plant diversity on fungal communities. However, plant richness and Shannon diversity were recorded 7 and 8 years before the experiment (in 2011 and 2012) (Fig. S5) and showed no differences between meadow and the cropping systems (organic and conventional treatments) but higher plant richness in organic compared to conventional treatment in 2012 (Fig. S5). Plant diversity has previously been positively linked to AMF diversity (Hiiesalu et al. 2014), indicating that high plant richness in organic treatment might partly explain why organic treatment did not differ from meadow in AMF richness whereas conventional treatment did. Higher plant richness in organic treatment is most probably a consequence of the lack of herbicide usage and is thus part of the management intensity effect. Organic and conventional treatment in this study already had a moderately diversified cropping system with 5-year rotation which included grass and crop mixtures (Salonen et al. 2023). However, decreasing management intensity by incorporating reduced tillage and increasing plant diversity by, for instance, cover-cropping, where noncommercial plants are grown together or after the main crop, could potentially further promote AMF richness in organic and conventional treatments (Thapa et al. 2021).
Arbuscular mycorrhizal fungal communities were affected by treatment and depth, but no treatment-specific taxa were found
We took a closer look at the AMF communities since the beneficial functions associated with AMF, such as induced nutrient uptake and protection against pathogens, can differ between AMF taxa (Sikes et al. 2010). We found AMF communities to be affected by treatment and depth but AMF taxa-specific differences between meadow, organic, and conventional treatments were not found. Based on patterns of fungal biomass allocation, AMF taxa can be grouped into rhizophilic guild, that have high biomass in roots and may protect host plants from pathogen colonization, edaphophilic guild, that have high extradical hyphae biomass and improve plant nutrient uptake (Weber et al. 2019), and ancestral guild, that produce low biomass both within and outside the root (Treseder et al. 2018, Phillips et al. 2019). In our study, the rhizophilic AMF guild was most pronounced, followed by the ancestral AMF guild. High proportion of rhizophilic AMF guild indicates an improved protection over plant pathogens. Edaphophilic guild was the least represented AMF guild in the studied soils, although the only edaphophilic genus, Diversispora, was found in all treatments. The abundance of many edaphophilic AMF taxa, but not Diversispora, has been linked to a higher C-N-ratio than what was present in our soils (Treseder et al. 2018, Fig. S4). Yet, the presence of a plant-nutrient-uptake improving AMF taxa such as Diversispora in organic and conventional treatment is an encouraging finding as it could benefit crop plants by scavenging large volume of soil, including deep soil, for nutrients.
In addition to depth, fungal trophic modes were affected by land use and agricultural management
In all soil treatments, the proportion of symbiotrophic fungi increased toward deeper soil layers, and in meadow, organic, and conventional treatment this was shown as an increase of AMF proportion in subsoil in comparison to topsoil. Since AMF benefit plant nutrient uptake (Smith and Smith 2011), and subsoil can harbour more than two-thirds of the nutrients in arable fields (Kautz et al. 2013), this subsoil association of AMF could indicate an important role of subsoil as a nutrient pool in the studied meadow, organic, and conventional treatments. Regarding forest treatment, our results support the previously proposed hypothesis that symbiotrophic mycorrhizal fungi in boreal forests are more competitive than saprotrophs in deeper layers where litter is more decomposed and C:N ratio is lower (Lindahl et al. 2007, van der Wal et al. 2013, Santalahti et al. 2016, Carteron et al. 2021), as both the highest symbiotrophic proportion and the lowest C:N ratios coincided in the same deep forest soil layers (30–40 cm and 40–80 cm) (Fig. 3; Fig. S4). Our results further suggest that the direct fungal interactions with plants, whether symbiotrophic or pathotrophic, are emphasized in deep forest soil (40–80 cm), where symbiotroph and pathotroph-saprotroph fungi represented the majority (75% and 18%) of fungal functional community and pure saprotrophs only a marginal (2%) (Fig. 3C). This indicates that the role of aboveground vegetation in shaping fungal communities in subsoil of boreal forest may be substantial.
Pathotrophic fungi were affected by treatment and depth. Out of all fungal functional groups, pathotrophic fungi correlated most strongly and negatively with depth (Fig. 3; Table S5), which may be explained by lower host interactions due to lower plant input and the typically lower richness of protist and soil animals in deeper soil layers (Du et al. 2022, Islam et al. 2022). Yet, contrary results were previously observed in a study with wheat-cropping system, where pathotrophic fungi were either unaffected or positively affected by depth (Schlatter et al. 2018). Among the pathotrophic fungi, plant pathogens were strongly influenced by treatment. We found organic and conventional treatments to increase plant pathogen proportion compared to forest and meadow (Table S6). Our results do not agree with a previous study where plant pathogen richness and proportion were shown to increase according to SOC (Du et al. 2022). Rather, our results are in line with the plant pathogen-inducing effect of arable soils over grasslands (French et al. 2017).
Saprotrophic fungi have been gaining attention as a potentially beneficial fungal group in agricultural soils contributing to nutrient cycling, soil fertility, plant pathogen suppression and SOC (Deacon et al. 2006, van der Wal et al. 2013, Ning et al. 2021, Hannula and Morriën 2022). We found the proportion of saprotrophic fungi to be more associated with low-intensity meadow treatment than with the cropping systems, organic and conventional treatments. Although the decomposing function of saprotrophic fungi can increase soil respiration and loss of carbon from soil in some cases (Newsham et al. 2018), a positive link between saprotroph biomass and SOC is frequently observed in agricultural soil (Six et al. 2006). In our study, higher SOC (Salonen et al. 2023) and higher saprotroph proportion coincided in meadow (Tables S6 and S8), further supporting the role of saprotrophs in SOC accrual.
Fe-ox were positively related to fungal communities down to 40–80 cm soil layer with strong correlation with arbuscular mycorrhizal fungal richness
Several soil properties contributed to fungal communities in meadow, organic and conventional treatment (Fig. 2; Tables 1 and 2). Fungal community differences (Bray–Curtis) were influenced by soil properties commonly observed in previous studies, C, N, DOC, C/N, P-tot, and pH (Francioli et al. 2016, Khan et al. 2016, Muneer et al. 2021, Rousk et al. 2010, 2011, Tedersoo et al. 2014, 2020, Zheng et al. 2019), as well as by root biomass, P-org, and Fe-ox, and less by P-inorg and P-H2O. Our results confirm the less commonly reported role of root biomass along the soil vertical profile in shaping fungal communities as well as the positive correlation of root biomass with fungal and AMF richness (Broeckling et al. 2008, Eisenhauer et al. 2017, López-Angulo et al. 2020). Previously, the role of P in shaping fungal communities has been emphasized, especially in agricultural soil (Francioli et al. 2016, He et al. 2016, Wu et al. 2022). Here, we consistently found P-org out of the different P forms (P-org, P-tot, P-inorg and P-H2O) to best explain variations in fungal community differences and fungal richness. Additionally, P-org explained fungal community differences better than any other soil property when the whole soil profile was considered. Total and available P has been shown to correlate negatively with fungal diversity (Wu et al. 2022), and in general, soil P is believed to negatively affect fungal richness (Tedersoo et al. 2014). In contrast, we did not find a negative link between fungal richness and any P form measured, and only a weak negative link with P-inorg and AMF richness was observed. Although the negative effects of soil P on AMF richness and abundance are well documented (Abbott et al. 1984, Camenzind et al. 2014, Chen et al. 2014, Jasper et al. 1979, Mosse 1973, Olsson et al. 1997), recently opposing effects of P in topsoil (negative) compared to subsoil (positive) were found (Luo et al. 2021), indicating the effects of P on AMF richness to vary within the soil vertical profile. This could explain why we did not find a negative link with most P forms and AMF richness when assessing the whole soil profile. However, our results do support the strong adverse role of different P forms on AMF proportion (Table 2) and suggest that AMF diversity and proportions may be differently affected by P in the soil vertical profile.
The strong role of Fe-ox in fungal communities is not commonly reported, making it a novel and interesting finding (Brandt et al. 2024, Jeewani et al. 2020). Fe-ox was the only soil property that associated with fungal communities consistently throughout the soil vertical profile (0–80 cm) and additionally correlated strongly with fungal and AMF richness (Tables 1, 2). This is supported by a previous study, where AMF were reported to preferentially associate with iron oxide surfaces in rhizosphere soil (Whitman et al. 2018). Fe-ox is important in soil aggregate formation and is associated with SOC (Jeewani et al. 2020, Pronk et al. 2011, Salonen et al. 2023), which can at least partly explain the role of Fe-ox as both soil aggregates and SOC are known to shape fungal communities (Fan and Wu 2021, Upton et al. 2019; Yang et al. 2019). Additionally, soil P availability is negatively affected by Fe-ox, as well as by Al-ox, which adsorb phosphate ions through ligand exchange reaction (Hingston et al. 1967). Thus, Fe-ox may have affected fungal communities by controlling the amount of available P. Further studies are needed to better understand the role and function of Fe-ox in shaping fungal, especially AMF, communities.
Meadow treatment had the most distinct soil properties among meadow–organic–conventional soil management intensity gradient
As such, pH influences fungal community structure (Hannula et al. 2021, Tedersoo et al. 2020) and it was overall higher in the cropping systems (organic and conventional treatment) compared to meadow (Table S8). Lower pH may also have led to higher Fe-ox content in meadow (Thompson et al. 2006). In addition to differences in pH and Fe-ox, meadow treatment was associated with higher C, N, C/N, DOC, Al-ox and root biomass in most soil layers, and higher P-org in the topsoil (0–10 cm) compared to the cropping systems, indicating their role in fungal community differences between meadow and conventional/organic treatment. Root biomass, C, N, P-org, and P-H2O were the only soil properties that significantly differed between organic and conventional treatments at least in some of the soil layers, indicating a link between these soil properties and variations in the fungal communities (Table S8). In deep soil, the role of root biomass may have been important as it was the only significantly different soil property between organic and conventional treatments in 30–40 cm and 40–80 cm soil layers. As Fe-ox differed significantly only between meadow and the cropping systems (organic and conventional), its role in the fungal community differences between organic and conventional treatments remains unclear. However, soil layer and treatment alone better explained the observed fungal community differences than all the measured soil properties together (PERMANOVA; R2 = 0.41 vs. R2 = 0.34). This indicates that soil management and depth may influence fungal communities beyond these commonly measured soil properties.
Conclusions
Our experimental set-up made it possible to study the long-term impacts of land use and soil management intensity on fungal communities. We showed that the effects of land use and soil management intensity on fungal communities persisted throughout the soil vertical profile down to 40–80 cm. In accordance with our hypothesis, the less intensively managed meadow was more associated with the potentially beneficial fungal groups than the more intensively managed organic and conventional cropping systems by having the highest AMF richness and saprotroph proportion. However, the management intensity differences between organically and conventionally managed soils were not reflected in significant differences in the potentially beneficial fungal groups. Organic and conventional treatments were distinguished by having the highest pathotroph richness and pathotroph and plant pathogen proportion, and forest by having the highest symbiotroph proportion. Similarly, as in forest but on a smaller scale, the mycorrhizal mode of requiring nutrients and energy became proportionally more important in deeper soil layers of meadow, organic, and conventional treatments indicating that subsoil nutrient reservoir could potentially be better utilized and the environmental impacts of farming reduced by optimizing agricultural soil management toward AMF favoring practices. We showed several fungal taxa to be proportionally more prominent in certain soil layers. Topsoil-associated taxa in meadow, organic, and conventional treatments, included the fungal classes Dothideomycetes, Sordariomycetes, and Tremellomycetes. Subsoil-associated taxa included the fungal classes Mortierellomycetes in all treatments and Leotiomycetes in meadow, organic, and conventional treatments. This study showed that the only soil property consistently significantly related to fungal communities throughout the soil vertical profile and with strong positive correlation with AMF richness was Fe-ox, which should be further studied. Additionally, this study indicated that sampling depth should be extended at least to 30 cm deep to better describe the diversity of AMF. Vertical profiles of agricultural soils that deploy more extensive regenerative agricultural practices, such as cover-cropping with deep-rooting plants and minimal-tillage, should in the future be explored for their microbial communities to better understand soil management effects in subsoils.
Acknowledgments
The authors wish to thank Riitta Lemola for the experimental setup, Matti Ylosmäki for sampling the soil, and Jenni Jääskeläinen, Eija Hagelberg, Stephanie Duranceau, Aku Pakarinen, Miia Collander, and Rashmi Shrestha for the help in soil sampling and laboratory analyses.
Author contributions
Laura Häkkinen (Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Visualization, Writing– original draft), Igor S. Pessi (Methodology, Supervision, Writing– review & editing), Anna-Reetta Salonen (Data curation, Formal analysis, Investigation, Methodology, Resources, Writing– review & editing), Oona Uhlgren (Data curation, Formal analysis, Investigation, Methodology, Resources), Helena Soinne (Methodology, Writing– review & editing), Jenni Hultman (Conceptualization, Methodology, Supervision, Writing– review & editing), and Jussi Heinonsalo (Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Validation, Writing– review & editing)
Conflict of interest
None declared.
Funding
This work was supported by the Strategic Research Council (SRC) at the Research council of Finland [327342, 352435, 352436, and 340169] and The Finnish Cultural Foundation [00230457].
