Liming amendment was effective for altering soil chemistry
Each of five treatment sites had 10 tonnes per hectare of crushed lime dropped from a helicopter. The amount of crushed lime that was measured to have fallen at each site varied between 3.3 and 12.2 tonnes per hectare (means of 3.3, 7.7, 5.43, 12.2 and 5.1 tonnes per hectare were measured for sites 1 to 5, respectively), reflecting heterogeneity in the application method. Soil chemistry measurements (pH, calcium (Ca2+), aluminum (Al3+), magnesium (Mg2+), and manganese (Mn)) were taken one year post liming treatment (10months before the first microbiome samples). No soil chemistry measurements were taken prior to liming treatment. Soil pH, and exchangeable (plant available) Ca2+, Mg2+, and Mn levels were significantly higher while Al3+ levels were significantly lower (T-test p \(\le\) 0.05) in treatment than control sites in the Upper and Lower Forest Floor horizons (Supp. Figure 1, Fig. 2). The magnitude of the differences between control and treated sites was larger in the Upper Forest Floor (Supp. Figure 1, Fig. 2). In the Upper B Horizon, only Mn levels were significantly higher in the treatment than control samples, with no other significant differences observed between control and treated sites. Taken together, these data demonstrate that the liming treatment was effective for: increasing soil pH, increasing calcium, magnesium, and manganese concentrations, and decreasing aluminum concentrations in the Upper and Lower Forest Floor samples (Supp. Table 1).
Prokaryotic but not fungal alpha diversity is increased by liming treatment
We compared the alpha diversity between our samples by grouping them by their treatment (control or treatment) and their soil horizon (upper forest floor, lower forest floor and upper B horizon) for both the bacterial/archaeal (16S rRNA gene) and fungal data (ITS2 gene). For bacteria/archaea, there is significantly (ANOVA p \(\le\) 0.05) higher Faith’s phylogenetic diversity in the treatment vs control (Fig. 3A) samples, and in the upper forest floor and upper B horizon in comparison to the lower forest floor (Fig. 3B). When examining differences between control and treatment in the different horizons, Faith’s phylogenetic diversity was higher in the treatment than control in both the upper and lower forest floor – but was only significant (ANOVA p \(\le\) 0.05) in the upper forest floor – while the trend was reversed in the upper B horizon, and this difference was not statistically significant (Fig. 3E). Similar trends were observed using other alpha diversity metrics: number of ASVs, Chao1 richness, Shannon diversity and Simpson’s diversity (Supp. Figure 2). In contrast there are no significant differences between control and treatment samples within the fungal community for any alpha diversity metrics, within any of the soil horizons (Fig. 3C, 3F and Supp. Figure 3). There is significantly higher (ANOVA p \(\le\) 0.05) richness (number of ASVs and Chao1 richness) in the upper forest floor than either the lower forest floor or upper B horizon samples, but alpha diversity indices that also incorporate abundance measurements or evenness showed no differences (Fig. 3D and Supp. Figure 3).
Shifts in microbial community composition (beta diversity) are seen with liming amendment in prokaryotic but not fungal communities
We compared the differences in microbial diversity with liming amendment as well as between the different soil horizons, using Weighted UniFrac distance for the prokaryotic community and Bray-Curtis dissimilarity for the fungal community, with PERMANOVA tests that also included site, sample within site, and season (Figs. 3G-3N). For the bacterial/archaeal community, most of the variation between samples was due to the soil horizon (R2 = 0.433, p = 0.001), although there were also significant differences with liming treatment (R2 = 0.034, p = 0.001) and season (R2 = 0.024, p = 0.003), as well as a number of significant interactions between these variables, including between treatment and soil horizon (R2 = 0.034, p = 0.008; Fig. 3G and Supp Tables 1 & 2). There is clear clustering of all upper B horizon samples from samples from the other horizons on both the first (36% variation) and second (25% variation) PCoA axes (Fig. 3G), and while there is clustering of the lower forest floor samples away from the upper forest floor samples, this difference is not as large. Within the upper forest floor, there is clearly distinct clustering of the treatment samples from the control (Fig. 3I). By contrast, there is no clear clustering of the treatment from control within either the lower forest floor (Fig. 3J) or upper B horizon (Fig. 3K). The betadispersion analysis indicates that while there are significant differences between groups (see above) there are also significant differences within the treatment (betadisp, ANOVA p ≤ 0.01) and soil horizon groups (betadisp, ANOVA p ≤ 0.005).
Within the fungal community there were no significant differences with any of the variables that we included in the PERMANOVA tests (Fig. 3H and Supp. Tables 1 & 2). Unlike for the bacterial/archaeal community, for the fungal community there was no single variable that accounted for a large proportion of the variation between samples (maximum R2 = 0.056 for soil horizon). Furthermore, we note that for the fungal ordination plots only a small amount of variation is shown on the first two axes (7% and 6% for PCoA1 and PcoA2, respectively), indicating that the variation between the samples is not well represented on two axes. We observed some clustering of the upper and lower forest floor samples away from the upper B horizon (Fig. 3H), indicating that – as for the bacteria/archaea – the upper B horizon is slightly distinct from both forest floor horizons. Additionally, we observe minor clustering of the treatment from the control samples within the upper forest floor (Fig. 3L) and lower forest floor (Fig. 3M) that is not seen in the upper B horizon (Fig. 3N). The dispersion of the treatment groups is significantly different between treatment groups (betadisp, ANOVA p ≤ 0.05) but not between soil horizon groups (betadisp, ANOVA p ≥ 0.05). Overall, the composition of the fungal communities is not impacted significantly by treatment and horizon.
Taxa associated with liming amendment of chronically acidified soils
We have examined the abundance of taxa using both relative abundance, as this is what most studies to date have used, and rCLR abundance, as this accounts for the compositionality of microbiome data43. For the relative abundances, higher values indicate that a larger proportion of the sequences within samples belong to a particular taxon. For the rCLR abundance, a zero value indicates that the abundance of a taxon is equal to the mean log2 abundance of all taxa, with positive or negative values indicating higher or lower abundances than the mean log2 abundance, respectively.
In the prokaryotic community, there were seven phyla that were present in all samples (Proteobacteria, Bacteroidota, Acidobacteria, Planctomycetota, Verrucomicrobiota, Actinobacteriota, and RCP2-54) and some of these were also high in relative abundance in a large number of samples: (i) Acidobacteriota dominate in all sample types (42–67% relative abundance); (ii) Proteobacteria are the second most abundant phyla in all sample types (15–29%); (iii) Bacteroidota are abundant in Upper (10%) and Lower (3.4%) Forest Floor Treatment samples; and (iv) Planctomycetota account for ~ 3-5.5% relative abundance in all samples, but are more abundant in Upper Forest Floor samples (Supp. Figure 4).When we look at lower taxonomic levels, there are a large number of ASVs that were not able to be classified; only 8,114 of the 13,684 16S rRNA gene ASVs had genus-level classifications, with a large number of these being made up of genera that are likely not well characterized, i.e., Candidatus classifications (431 ASVs), or classifications without Latin scientific names (3,454 ASVs).
The fungal community was dominated by the Basidiomycota (78–90% relative abundance) and the Ascomycota to a lesser extent (8.7–20%); these were the only two phyla that were present in all samples (Supp. Figure 5) and all other phyla were present at < 2% relative abundance on average. The Mortierellomycota were also highly prevalent (present in 95% of Lower Forest Floor Control samples, and 100% of all other sample groups) but were only present at a maximum relative abundance of 1.7%. In contrast to the bacterial/archaeal community, where almost all taxa that were abundant were also highly prevalent, in the fungal community there were very few taxa with genus-level classifications that were present in > 90% of samples within at least one sample grouping; these belonged to the Agaricomycetes (Galerina, Sebacina, Amanita, Russula and Cortinarius), Mortierellomycetes (Mortierella) and Dothideomycetes (Cenococcum). For the fungal ASVs there was also a large proportion without a genus-level classification (only 1,210 of 1,785 ASVs had a genus-level classification).
Identifying differentially abundant fungal genera between control and treatment samples
To identify which genera might be driving the differences in community composition seen in the ordination plots (Fig. 3), we carried out differential abundance testing between control and treatment samples separately for each soil horizon. These tests were carried out using ALDEx2, ANCOM-II and MaAsLin2 and we consider a genus to be differentially abundant if two of these three tests identify it (see Methods section for further details). Overall, seven genera were identified as significantly differentially abundant, predominantly within the upper forest floor horizon, however, most of these were found in very low abundances. If we focus on only the 25 most abundant genera by relative or rCLR abundance (44 total), only three genera were significantly differentially abundant within these (Fig. 4): Apodus (more abundant in treatment samples in the upper forest floor horizon), Galerina (more abundant in control samples in the upper forest floor horizon) and Tylospora (more abundant in control samples in the upper B horizon; these trends were consistent between both relative and rCLR abundance).
Bacterial/archaeal taxa associated with lime amendment.
For the bacterial/archaeal genera we again carried out differential abundance testing between control and treatment samples within each soil horizon using ALDEx2, ANCOM-II and MaAsLin2 and we consider a genus to be differentially abundant if two of these three tests identify it. Overall, 230 genera were identified as being significantly differentially abundant between control and treatment samples in at least one of the three horizons; 185 in the upper forest floor, 13 in the lower forest floor and 8 in the upper B horizon. We focused on those that were also within the most abundant genera, taking the top 25 most abundant genera by relative or rCLR abundance (46 total; Fig. 5). Reflecting the magnitude of alpha and beta diversity shifts in the horizons, 24 or three genera (of the 46 total) were significantly differentially abundant between control and treatment samples in the upper or lower forest floor samples, respectively. Three of these taxa (all unclassified at the genus level; Commomonadaceae, Bacteroidia and Microscillaceae) were identified in both upper and lower forest floor samples, and are all more abundant in treatment than control samples. One genera was significantly differentially abundant in the upper B horizon; Acidipila (Acidobacteriae), slightly more abundant in treatment than control samples, 0.41 vs 0.18% and − 1.66 vs -2.05 for rCLR abundance, respectively. Five genera belonging to the Acidobacteriae class were identified as differentially abundant by at least 2 tests and are more abundant in control vs treatment samples. In contrast, the five genera identified by at least two tests belonging to the Bacteroidia class are more abundant in treatment vs control. Similar trends are not seen within any other classes that have multiple genera identified as differentially abundant.