Fungal community and functional responses to soil warming are greater than for soil nitrogen enrichment

Forest soils are an important and large carbon (C) sink (Pan et al., 2011). They harbor hyperdiverse fungal communities that act as key decomposers (Schneider et al., 2012; Peay et al., 2016) and, in turn, are regulators of forest soil C storage (Averill and Hawkes, 2016; Frey, 2019; Lindahl et al., 2021). Increasing evidence shows that fungal communities are sensitive to environmental changes, including anthropogenic nitrogen (N) deposition (Lilleskov et al., 2011; Morrison et al., 2016; Van der Linde et al., 2018; Moore et al., 2021), climate warming (Treseder et al., 2016; Fernandez et al., 2017; Morrison et al., 2019; Cao et al., 2020), nonnative invasive species (Lekberg et al., 2013; Gibbons et al., 2017; Anthony et al., 2019), and myriad other global changes (see review by Zhou et al., 2020). How fungi respond to global changes at the local level can have cascading effects throughout forests, influencing tree growth and mortality and soil carbon storage. However, and until recently, it has been technically challenging to trace how fungal community compositional changes are linked to explicit fungal functional shifts which impact processes such as soil C storage.

Ectomycorrhizal fungi (EMF) are tree symbionts that are especially important mediators of organic N acquisition from soil organic matter (SOM) and can strongly affect overall SOM chemistry (Read and Perez-Moreno, 2003; Lindahl and Tunlid, 2015). EMF relative abundances in forest soil are positively correlated with proteolytic enzyme activities (Averill and Hawkes, 2016; Morrison et al., 2019), soil C:N ratio (Anthony et al., 2017), and soil N immobilization (Franklin et al., 2014; Clemmensen et al., 2015). There is also evidence for a narrow set of EMF species to drive SOM cycling within forests (Clemmensen et al., 2015; Heinonsalo et al., 2015; Lindahl et al., 2021). For example, the relative abundance of Cortinarus acutus, an EMF species capable of producing manganese peroxidases, was linked to a 33% decrease in C stocks in the organic soil horizon of a Swedish boreal forest (Lindahl et al., 2021). Thus, predicting how forest SOM will vary with future global changes can be improved with concomitant knowledge of the fungal community, including ectomycorrhizal species.

Climate warming and atmospheric N deposition are particularly important environmental factors that can influence soil fungi, SOM cycling, and soil C storage. Warming often decreases fungal biomass (DeAngelis et al., 2015; Morrison et al., 2019) and has been shown to shift fungal community structure (Morrison et al., 2019; Cao et al., 2020) and the relative abundances of key fungal taxa (Fernandez et al., 2017; Mucha et al., 2018). Warming can also alter SOM chemistry by accelerating decomposition, resulting in lower soil C storage (Pisani et al., 2015; Melillo et al., 2017). Simulated N deposition also reduces fungal biomass at high N addition levels (Moore et al., 2021), shifts fungal community structure (Lilleskov et al., 2011; Morrison et al., 2016), and alters SOM chemistry by inhibiting decomposition. Most previous studies have focused on one or the other of these global change drivers, yet these are happening simultaneously in the real world. Because fungal communities often exhibit high degrees of functional redundancy (Rineau and Courty, 2011; Banerjee et al., 2016), it is also uncertain whether shifts in fungal diversity and community composition in response to warming and N enrichment ultimately translate into functional changes which impact SOM decomposition.

In this study, we measured how the diversity and composition of dominant fungal functional groups (i.e., EMF, arbuscular mycorrhizal fungi [AMF], saprotrophs, pathogens) responded to soil warming and N addition, alone and in combination. We also measured fungal functional genes using a gene capture approach which allowed us to target specific aspects of fungal decomposition which may be sensitive to environmental change and linked to SOM chemistry at the molecular level. Warming at our study site has enhanced soil respiration by 30%–50% (Contosta et al., 2011) and reduced C storage by 35% (Anthony et al., 2020). Nitrogen additions have had ephemeral effects on soil respiration and no significant effect on soil C storage. Warming in combination with nitrogen has had the greatest positive effect on soil respiration and intermediate effects on soil C storage (unpublished data). As the primary decomposers in this ecosystem, we expected shifts in the fungal community to mirror changes in soil respiration and C storage. Specifically, we hypothesized that warming and warming × N would have larger effects on fungal communities and decomposition genes compared with N addition alone. All of the treatments have shifted SOM chemistry to varying degrees (Pisani et al., 2015), with a general increase in the decomposition of most SOM compounds in the heated plots and a suppression of SOM decomposition in the N addition plots (vandenEnden et al., 2021). Because warming in combination with N has shifted SOM chemistry to be more similar to the heated compared to the N addition plots (vandenEnden et al., 2021), we expected fungal functional gene changes in the heating × N plots to be more similar to heated plots and distinct from the N addition plots. Earlier studies on fungi at this site were performed on a subset of the experimental plots (i.e., control vs. heated only plots) during the initial study phase (Pec et al., 2021) and in the context of an encroaching nonnative plant (Wheeler et al., 2017; Anthony et al., 2020), but the sensitivities of soil fungi to long-term warming, N additions, and their combination independent of plant invasion have not been explicitly investigated.

The Soil Warming × Nitrogen Addition experiment at the Harvard Forest Long-Term Ecological Research site (42°50′ N, 72°18′ W) was established in 2006 and is described in detail in Contosta et al. (2011). It is located in a mixed deciduous stand with Quercus rubra, Quercus velutina, Acer rubrum, Acer pensylvanicum, Fagus grandifolia, and Betula papyrifera canopy trees. Mean annual temperature is 8.3ºC and mean annual precipitation is 1,247 mm (Boose, 2021). The full experiment includes 3 × 3 m plots replicated 6 times for each of 4 treatments: control, heated, N addition, and heated × N addition. Heated plots are elevated 5ºC above ambient soil temperatures using buried heating cables installed at 10 cm depth and spaced 20 cm apart; N addition plots receive equal monthly doses of aqueous ammonium nitrate (0.83 g N m^–2) for a total of 5 g N m^–2 yr^–1 throughout the growing season (May–October).

We collected soil samples in July 2016 from 5 of the 6 replicate plots from each treatment. Five rather than 6 replicates were sampled because 1 heated plot is no longer functional due to an electrical issue, and we chose to have a balanced design. Within each plot, we collected a 10 × 10 cm block from the organic horizon to the depth of the mineral soil, followed by collection of a cylindrical core of mineral soil (5 cm width × 10 cm depth). Two paired organic and mineral soil samples were collected within each plot and homogenized into composites by depth increment (i.e., organic horizon and 0–10 cm mineral soil). Samples were stored on ice in the field and then at 4°C within 12 h of sampling. Within 24 h, soil was passed through a 4 mm sieve, and a subsample for molecular analysis was placed at –80°C. Additionally, a subsample for phospholipid and neutral lipid fatty acid analyses to measure fungal and AMF biomass, respectively, was frozen at –20°C and then freeze dried. All remaining soil was dried at 105°C for 48 h for analysis of soil pH (10 g soil: 20 mL deionized water) and total soil organic C and N using dry combustion (Perkin Elmer 2400 Series II CHN, Waltham, MA). A summary of the edaphic responses to the treatments can be found in Anthony et al. (2020).

We characterized fungal diversity, community composition, the relative abundance of key functional groups, and fungal functional genes related to decomposition. Genomic DNA was extracted from 250 mg of soil (Qiagen PowerSoil Kit; Qiagen, Hilden, Germany). Total fungi, including saprotrophs, EMF, and pathogenic species were analyzed using ITS2 DNA metabarcoding with the primer pair fITS7-ITS4 (White et al., 1990; Ihrmark et al., 2012; BioProject: PRJNA522440), and AMF were analyzed using 18S DNA metabarcoding with the primer pair NS31-AML2 (Simon et al., 1992; Lee et al., 2008; PRJNA522442). A description of the polymerase chain reaction (PCR) and thermocycler conditions is found in Anthony et al. (2020). Fungal functional genes were targeted and enriched using biotinylated RNA probe capture. A set of 20,005 probes were designed by Arbor Biosciences from 2,322 fungal gene sequences (a list of the gene sequences can be found in Moore et al., 2021). Each probe was 100 bp and varied with 3× tilling density. Probes were hybridized with sheared template DNA (350 bp) on a thermocycler using the manufacturer protocol and the MyBaits kit (Arbor BioSciences, Ann Arbor, MI, USA). Hybridized DNA was isolated and amplified using PCR. Amplified libraries were prepared for sequencing using the PCR-free TruSeq Prep Kit (Illumina, San Diego, CA, USA) and were sequenced on an Illumina Next-Seq 150 platform (PRJNA63326). A thorough description including additional details on the probe designs and thermocycler conditions is in Moore et al. (2021).

All sequences (ITS, 18S, functional genes) were passed through quality control measures, removing poor quality reads (Phred scores <2 at a 20 bp sliding window) using Trimmomatic (Bolger et al., 2014). Forward and reverse reads were then merged at a 20 bp overlap with a 5% mismatch (ITS2) and at a 10 bp overlap allowing 10% mismatch (18S) using the join_paired_ends.py function in QIIME (Caporaso et al., 2010). The entire ITS2 region was then excised from flanking 5.8S and 28S regions using ITSx (Bengtsson-Palme et al., 2013). Quality control sequences were then clustered into operational taxonomic units (OTUs) using USEARCH v8 (Edgar, 2010) at a 97% sequence similarity, removing singletons and chimeras with the cluster_otus function. ITS2 OTUs were assigned taxonomy using the UNITE database v7 (Abarenkov et al., 2010) and the assign_taxonomy.py function in QIIME; 18S OTUs were blasted against the MARJAAM database (Öpik et al., 2010) and then the National Center for Biotechnology Information (NCBI) nt database for nonassigned sequences. Nonfungal OTUs were removed from all data sets and non-Glomeromycotina OTUs were removed from the AMF data set. The taxonomy of ITS2 OTUs of statistical significance based on similarity percentage analysis (see below) were manually curated using BLAST searches against the UNITE database, accepting the best hit (lowest E value) if sequence identity was >98%. ITS2 reads were also assigned potential functional guild annotations using FUNGuild (Nguyen et al., 2016), accepting all “probable” or higher level assignments. For ectomycorrhizal and pathotrophic fungi, a taxon was considered ectomycorrhizal or pathogenic even if it was also identified as possessing >1 trophic mode (e.g., an ectomycorrhizal/ericoid assignment was considered ectomycorrhizal).

Quality-controlled functional genes were aligned to a custom database made from the NCBI template DNA sequences originally used to design the probes using bwa v7 (Li, 2013). The database was made using the index algorithm, and forward and reverse reads were aligned using the bwa mem algorithm. We only aligned forward sequences since forward and reverse sequences did not pair well, and we removed alignments with <90% sequence similarity using SAMtools. Sequences were then sorted by genes encoding for hydrolytic enzymes (betaglucosidase, cellobiohydrolase, and cellobioside dehydrogenase) or oxidative enzymes (lignin peroxidase, manganese peroxidase, laccase, and laccase-like multicopper oxidases).

Molecular components of SOM were analyzed in soil extracts using gas chromatography-mass spectrometry (GC-MS) with an Agilent 7890B GC with a 5977B MS with electron impact ionization. See vandenEnden et al. (2021) for a complete description of the methods. In short, we quantified n-alkanes, n-alkanoic acids, sugars, cutin-, suberin-, and lignin-derived (vaniyllyl, syringyl, and cinnamyl phenols) compounds. Compounds were identified using a library of MS spectra (Wiley registry v9), National Institute of Standards and Technology (2008), and a custom mass spectral library. Quantification was based on standards specific to each compound class and normalized by soil C content (mg g^–1 C).

All statistical analyses were conducted in R (v3.6.1; R Development Core Team, 2019) and significance was set to a P value ≤ 0.05. We used analysis of variance (ANOVA) to test the effects of warming and N additions plus their interaction on univariate response variables using the base R aov function. We used the Anova function in the car package (Fox et al., 2007) to calculate type III (sequential) ANOVA tables. Multiple comparisons were evaluated by computing estimated marginal means using the emmeans function in the lsmeans package (Lenth, 2018). All models were inspected to meet the assumptions of ANOVA based on the distribution of residuals and QQnorm plots. All molecular data sets were normalized to the lowest sequencing depth using the rarefy function in vegan (Oksanen et al., 2013). We computed species richness and diversity (Shannon-Index) using the specnumber and diversity functions in vegan, respectively. We calculated Pielou’s evenness index J’ as Shannon/log(richness). We analyzed differences in community composition using permutational ANOVA (PERMANOVA) and the adonis function in vegan and visualized differences by treatments using non-metric multidimensional scaling (NMDS) and the metaMDS function in vegan. OTUs contributing most to dissimilarity between control and treatment plots were identified using similarity percentage analysis and the simper function in vegan. All multivariate analyses were based on fungal relative abundances computed as Bray–Curtis dissimilarities using the vegdist function in vegan.

Fungal co-occurrence was computed using probabilistic models and the cooccur function in the cooccur package (Griffith et al., 2016). We separated the data set based on significant differences in community composition to avoid spurious negative co-occurrences due to differences in environmental conditions (Blanchet et al., 2020). We then removed any species with fewer than 2 occurrences in the data sets. We used a combinatorics approach to compute co-occurrence (Veech, 2013) by setting prob = “comb.” At the plot level, we summarized every nonrandom, positive and negative pairwise co-occurrence that was observed. We then calculated our 2 response variables, the mean probability of positive and negative species co-occurrences at the plot level in order to compare overall co-occurrence probabilities across the treatments.

We tested the correlation between fungal community characteristics and functional genes with SOM compound chemistry and soil organic C using a regression framework. We first examined overall differences in SOM compounds using PERMANOVA and NMDS as described above. We then correlated fungal predictor variables, including community composition, represented as NMDS1 and NMDS2 for total fungi, EMF, saprotrophs, and pathotrophs, richness and Shannon diversity (all groups listed above), co-occurrence (negative and positive), and the relative abundances of fungal functional genes. We used the envfit function in vegan to examine correlations between fungal characteristics and SOM composition and linear regression to test the relationship between SOM compounds and every fungal predictor variable. Separate linear regression models were then run for every combination and inspected for normality. We used the log concentration of SOM compounds where nonlinear relationships existed.

Fungal biomass was not affected by the treatments (Table S1), but heating significantly altered fungal diversity and community composition while N additions only significantly affected fungal diversity. Fungal community composition in the organic horizon was distinct in ambient compared to heated plots (F_1,19 = 1.45, P = 0.03; Table 1, Figure 1), but this effect was not detected in the mineral horizon (Table S2). A similar response to heating in the organic horizon was observed for saprotrophic fungi (F_1,19 = 1.85, P = 0.01), EMF (F_1,19 = 1.42, P = 0.05), and AMF (F_1,19 = 2.32, P = 0.03; Table 1, Figure 1). In contrast to community composition, both heating and N additions reduced fungal (F_1,19 = 4.97, P = 0.04) and EMF Shannon diversity (F_1,19 = 12.16, P = 0.003; Figure 2, Table S3) and community evenness, measured as Pielou’s evenness index J’ (Table S4). Saprotrophic diversity was not sensitive to the treatments. AMF diversity (F_1,19 = 7.77, P = 0.009) and community evenness (F_1,19 = 9.5, P = 0.004) were also reduced by heating but not by N additions. Overall, fungal communities in the control plots were significantly more even than the treatment plots, and evenness contributed more to differences in diversity than richness. Changes in community composition and diversity across the treatments were observed in the organic horizon and to a lesser extent or not at all in mineral soil; thus, all subsequent analyses report on communities found in the organic horizon where we see the greatest impact on soil C pools.

EMF were the dominant fungal guild (approximately 70% relative abundance; Figure S1), and there were no treatment-level differences in EMF relative abundances as a group, nor that of any other fungal guild (i.e., saprotrophs, pathotrophs; Table S5). Rather, there were particular ectomycorrhizal, and to a lesser extent, saprotrophic OTUs that drove dissimilarity in fungal community composition between treatment and control plots (Table S6). Heated plots were hyperdominated by ectomycorrhizal Russula laurocerasi, while N addition plots were hyperdominated by Russula subsulphurea compared to control plots (i.e., these OTUs were at an order of magnitude higher relative abundances in treatment versus control plots; Figure 3; Table S7). In the 2-factor heated × N plots, ectomycorrhizal Boletus rubropunctus had 6 times greater relative abundance compared to control plots. Representing the core sensitive taxa, ectomycorrhizal Amanita fulva, Cortinarius sp. 171, and Cortinarius sp. 1565 and the nonmycorrhizal Eurotiales sp. 18, Herpotrichiellaceae sp. 143, and Mortierella pulchella were all negatively affected by the treatments.

Fungal community co-occurrence patterns were altered by heating but not by N additions. In both the heated and the heated × N addition plots, there was lower negative co-occurrence among fungal community members compared to unheated plots (F_1,19 = 4.4, P = 0.05; Figure 4A). Conversely, there was higher positive co-occurrence among fungi in the heated and heated × N addition plots (F_1,19 = 6.9, P = 0.02; Figure 4B). Taxa that most negatively co-occurred with other species in the heated plots included saprotrophic, EMF, and pathotrophic fungal species. There were no changes in the co-occurrence of AMF communities.

Overall fungal functional gene composition was not affected by the treatments, but the relative abundances of individual decomposition genes and groups of genes were altered by heating and N additions (Table S8). Fungal genes encoding for manganese peroxidases were at reduced proportions in the heated plots (Figure 5A), particularly in the single factor heated (i.e., heating only) treatment. In contrast, the total relative abundance of hydrolytic enzyme encoding genes was higher in the heated plots (Figure 5B). Although there was no change in the relative abundance of hydrolytic enzyme encoding genes in the N addition plots, cellobiohydrolase encoding genes were at reduced proportions (Figure S2).

Finally, we examined which fungal community attributes were associated with changes in SOM compounds and soil C storage in the treatment plots. Overall, SOM chemistry was altered by heating (PERMANOVA: F_1,19 = 4.5, P = 0.005) but not by N additions (F_1,19 = 1.3, P = 0.28), and the heated plots were enriched in suberin- and cutin-derived compounds compared to unheated plots (Figure 6A). The relative abundance of fungal hydrolytic enzyme encoding genes was positively associated with both suberin- and cutin-derived compounds. In contrast, lignin-derived compounds, short- and long-chain alkanoic acids, n-alkanes, and total sugars were associated with the unheated plots. Neither fungal community composition nor diversity (total or that of any group) were correlated with total soil C or any particular group of SOM compounds; however, fungal negative co-occurrence was strongly and negatively correlated with suberin- and cutin-derived compounds (Figure S3). Fungal negative co-occurrence was also strongly and positively correlated with overall soil organic C contents (Figure 6B).

Few studies have considered how soil fungi in forests will be affected by future warming conditions, atmospheric N deposition, and their combination (Rillig et al., 2019; Zhou et al., 2020). But as the main decomposers of SOM in forests (Schneider et al., 2012), fungal responses to warming and N additions may be linked to soil C storage in the face of climate change. Here, we show that fungal communities are less diverse (i.e., lower evenness), shift in composition, and have lower manganese peroxidase and higher hydrolytic enzyme production potentials in response to heating in comparison to N additions. This indicates that fungal potentials to decompose cellulolytic versus ligninolytic substrates responded in opposite directions to heating. The novel gene capture and enrichment technique we applied to soil systems made it uniquely possible for us to measure these specific fungal-mediated decomposition pathways. Changes in fungal functional genes and co-occurrence were also associated with distinct SOM compound concentrations and soil C content. Below, we discuss these results in greater detail, along with the potential implications for soil C storage.

Fungal diversity and community composition shifted in response to heating. Nitrogen additions had a small and insignificant effect on fungal community composition and a milder but significant effect on fungal diversity. Based on changes in soil respiration and C storage in the treatment plots, we expected fungal communities to be most affected by heating and heating × N followed by N additions alone, and our results generally support this prediction. Our findings are also consistent with previously described warming effects on fungal community composition at the Prospect Hill Soil Warming experiment also located at the Harvard Forest (Morrison et al., 2019; Pec et al., 2021). Greater warming versus N addition effects on fungal community composition have also been observed in a subtropical forest (Cao et al., 2020). In Cao et al. (2020) and our study, larger heating versus N addition effects may be due to the quantity and duration of added N. Because higher N additions (15 g N m^–2 yr^–1) shifted fungal community composition after more than 20 years of fertilization at a neighboring experiment within the Harvard Forest (Morrison et al., 2016), insignificant trends observed here are likely due to the smaller quantity and shorter duration of N additions (Moore et al., 2021). Heating was therefore a stronger factor than N additions affecting fungal community composition in our study, but the magnitude of this effect differed across fungal functional groups and with simultaneous heating and N additions. Notably, the only significant interaction between heating and N additions was observed for saprotrophic fungal community composition which shifted most in the heated × N addition compared to control plots. Interestingly, the heated × N addition plots also have the highest soil respiration rates. One possibility is that soil respiration may be linked to changes in saprotrophic fungal community composition to a greater extent than other functional groups.

There were no treatment effects on fungal biomass or the relative abundances of fungal trophic guilds. Fungal communities from all plots were dominated by mycorrhizal fungi followed by saprotrophs. Although at the same relative abundances, mycorrhizal diversity decreased with the treatment compared to control plots. In particular, community evenness declined while species richness remained unchanged. The ectomycorrhizal Russula laurocerasi and Russula subsulphurea were at more than an order of magnitude higher relative abundance in the heated and N addition plots compared to control plots, respectively, and were thus major contributors to observed declines in fungal evenness. A dominant AMF Glomus sp. that drove observed differences in community composition between the heated and control plots (Table S10) had more than 2 times greater relative abundance in the heated (45% relative abundance) versus control plots (18%; Figure S4). These observations of overdominance are consistent with other global change studies at the Harvard Forest showing that N and manganese additions select for hyperdominance by Russula vinaceae (Morrison et al., 2016) and Coccinonectria rusci (Whalen et al., 2018), respectively. Our results lend support to the idea that soil warming and N addition encourage hyperdominance by a handful of fungal species.

The functional implications of hyperdominance in mycorrhizal fungal communities remain largely unknown. One possibility is that diversity stabilizes communities to additional environmental changes (Elton, 1958; Hillebrand et al., 2008). For example, using synthetic communities, denitrification by uneven microbial communities shifted more in response to salinity stress than denitrification by even communities (Wittebolle et al., 2009). Soil fungi in the heated plots where our work was conducted were also more sensitive to the additional pressures of a nonnative, invasive plant, Alliaria petiolata (Anthony et al., 2020). Warming-induced reductions in fungal community evenness may reduce community resistance to additional environmental changes and may have myriad additional implications for forest functioning in a changing world.

Fungal taxa exhibited reduced negative and enhanced positive co-occurrence in the heated plots. This suggests that potential species–species interactions shifted with heating, favoring coexistence over competition. Saprotrophic fungi are widely known to compete for space and nutrients, and protagonists even actively preclude colonization by antagonists when decomposing wood (Heilmann-Clausen and Boddy, 2005; Boddy and Hiscox, 2017); but these interactions are challenging to observe and quantify at the microscopic level in soil. This is why an increasing number of studies use co-occurrence analysis to identify potential species–species interactions in DNA-based surveys (Berry and Widder, 2014; Ovaskainen et al., 2017; Liu et al., 2019). The challenge with this approach is that it is typically unclear whether co-occurrence patterns are due to species interactions (Gotelli and McCabe, 2002), stochastic drift processes (Ulrich, 2004), or shared environmental preferences (Ovaskainen et al., 2017). Because we deliberately made separate co-occurrence models for the heated and unheated plots, spurious direct effects due to differences in fungal community composition, soil C contents, and SOM composition were avoided. While we cannot fully tease these scenarios apart in our study, these co-occurrence patterns establish the possibility that fungal competition was altered by heating.

An immediate question is therefore why fungi might compete less in response to a decade of heating. One possible explanation is a C limitation hypothesis. Because soil warming progressively reduces soil C content, fungi may have become more cooperative (i.e., sharing resources) over time with warming. Simulation studies suggest that low C availability promotes microbial cooperation through the sharing and more efficient use of decomposition products (Ebrahimi et al., 2019). Experimentally, low C content can slow microbial growth (Sawada et al., 2008; Reischke et al., 2014) and select for species that produce fewer extracellular enzymes (Malik et al., 2018), creating the conditions for lower direct and indirect competition. Saprotrophic fungi with high competitive abilities also exhibit lower tolerance to temperature and moisture changes (Maynard et al., 2019), which supports the idea that heating might select for species with lower competitive abilities. The stress gradient hypothesis also states that facilitative versus competitive interactions should increase in frequency when physical conditions are suboptimal (Maestre et al., 2009), especially under conditions of energy and nutrient limitations, as observed for soil C and most SOM constituents in the heated versus unheated plots (Figure 6A). Thus, there is a theoretical basis for C limitations to reduce fungal competition and promote facilitative interactions consistent with our observations based on co-occurrence analysis.

Soil warming typically increases soil respiration via accelerated decomposition (Melillo et al., 2017; Romero-Olivares et al., 2017). While this positive effect of warming can diminish over time, at the Soil Warming × Nitrogen Addition Study where this work was conducted, soil respiration remained elevated compared to control plots after 10 years of continuous heating (unpublished data). As the primary decomposers in forests, changes in fungal-mediated decomposition may drive long-term soil C losses due to heating. However, it has been technically challenging to study fungal functions since currently available approaches do not capture high resolution fungal functional genes (e.g., metagenomics that mostly capture bacterial genes) or differentiate between fungal and bacterial decomposing enzymes (e.g., potential extracellular enzyme activities). Here, we used a novel gene capture technique (see Anthony et al., 2020; Moore et al., 2021) to specifically study how changes in fungal decomposition genes shifted with heating and N additions. While we expected fungal functional gene shifts to be the greatest and similar between the heated and heated × N compared to N addition only plots, we found only partial support for this with respect to manganese peroxidase and hydrolytic enzyme encoding genes.

Heating reduced fungal manganese peroxidase gene proportions while N additions had no effect. Because manganese peroxidases are often the most common lignin-decay enzymes produced by fungi (Entwistle et al., 2018), they can have large impacts on soil C storage. Earlier work at this experiment demonstrated that heating reduces lignin-derived compounds (Pisani et al., 2015). Thus, one possibility is that lower availability of ligninolytic substrates reduced fungal production of manganese peroxidases; however, there was no correlation between lignin-derived compounds and manganese peroxidase gene proportions. It is also possible that this change is due to reductions in soil manganese availability (Whalen et al., 2018). A recent meta-analysis found that soil warming in forests increases woody plant manganese concentrations by approximately 40% (Yuan et al., 2018), which could dramatically reduce soil manganese levels by immobilizing it in woody biomass. Overall, soil fungi in the heated plots have a reduced capacity to synthesize manganese peroxidases, an indicator of lower ligninolytic decomposition potential by fungi. Interestingly, long-term warming can select for bacteria more adapted to degrading lignin (Pold et al., 2015), and thus there may be a shift in the relative importance of fungi versus bacteria in lignin decomposition under warming conditions.

Conversely, heating increased the relative abundances of fungal hydrolytic enzyme encoding genes. This is consistent with a recent meta-analysis showing that warming generally increases hydrolytic enzyme activities and that this increase in hydrolytic enzymes is strongly correlated with fungal biomass (Meng et al., 2020). These enzymes act on cellulolytic substrates which may be reduced in the heated plots at our study site. Recently, vandenEnden et al. (2021) found lower O-alkyl C compounds in heated soil at our experiment using solid-state ¹³C NMR spectra. This is an indicator of lower cellulolytic compound content and may be linked to the increased fungal hydrolytic enzyme gene potentials identified here. Interestingly, shifts in hydrolytic enzyme gene proportions also mirror the trajectory of soil respiration rates at this experiment, consistent with our initial expectations. In particular, both respiration and hydrolytic enzyme gene proportions are elevated in the heated plots compared to control plots. In the N addition plots, neither respiration nor hydrolytic enzyme gene proportions are different from the control plots. In the heated × N addition plots, both respiration and hydrolytic enzyme potentials are elevated to an even greater extent than in the heated plots (unpublished data). Thus, shifts in fungal hydrolytic enzyme encoding genes and changes in soil respiration may be linked across the treatment plots. While we found a small decrease in cellobiohydrolase gene proportions with N additions, there was no overall change in total hydrolytic enzyme gene proportions, suggesting that this did not constrain the liberation of cellobiose from cellulose polymers (Rani Singhania, 2011). Regional studies over gradients of atmospheric N deposition show that hydrolytic enzyme encoding genes increase (Moore et al., 2021), but we did not verify this pattern with our experimental N addition. Our study shows that warming specifically increases fungal hydrolytic enzyme encoding gene proportions and that fungi may be key mediators of increased hydrolytic enzyme activities in heated soils.

The fate of soil C in a warmer world has important feedbacks to climate change. Since fungi are the dominant decomposers in forest ecosystems (Schneider et al., 2012), it was surprising that SOM compounds and C content were not strongly correlated with fungal diversity, community composition, nor that of any individual functional guild. However, fungal hydrolytic enzyme gene proportions were negatively correlated with soil organic C content and positively correlated with cutin- and suberin-derived compounds. Higher proportions of fungal hydrolytic enzyme encoding genes in the heated plots may indicate greater hydrolysis of SOM with implications for soil C storage. Thus, one possibility is that fungal functional shifts, more so than compositional changes, help to explain reduced soil C storage in chronically heated soil.

Striking above and beyond fungal functional gene changes was the positive correlation between fungal negative co-occurrence and soil C content. As discussed above, negative co-occurrence may be a reflection of fungal–fungal competition which was reduced in the heated plots. Increasing evidence suggests that specific forms of microbial competition may trade-off with other functions such as growth and the acquisition of particular resources (Maynard et al., 2017; Malik et al., 2019; Maynard et al., 2019). The production of antimicrobial toxins by fungi can constrain other metabolic functions (Czárán et al., 2002). For example, EMF species that are highly competitive for space on tree roots (e.g., Thelephora terrestris) have been linked to lower hydrolytic enzyme activities, whereas less dominant species (e.g., Suillus pungens) are associated with higher hydrolytic enzyme activities (Moeller and Peay, 2016). Thus, the decomposition of hydrolysable SOM constituents may be constrained by competing metabolic demands associated with competition. If fungi compete less in heated soil, then they may be free to invest in other processes, such as the hydrolysis of particular SOM constituents, which can reduce soil C content. Higher hydrolytic enzyme encoding gene proportions in the heated plots lends support to this idea. A presumed exception to a competition-decomposition trade-off is indirect competition for resources liberated via oxidative decomposition.

Of course, the processes underpinning soil C loss with long-term warming are more complex than can be understood from a single sampling campaign (Melillo et al., 2017). In this study, we are likely only capturing a unique experimental period characterized by fungal functional and species coexistence shifts linked to a particular phase of soil C loss. It is also important to consider what the directionality and direct versus indirect connections may be between fungal shifts and soil carbon storage, which cannot be disentangled here.

One unanswered question worth investigating in the future is whether C limitation reduces fungal competition or if soil C storage decreases when fungi compete less and decompose more. While it may be challenging to answer these questions in the field, this could be addressed experimentally. For example, using competition mesocosms involving negatively co-occurring taxa (sensu Maynard et al., 2017), it would be possible to measure decomposition in the presence or absence of competition and under manipulated C concentrations. Despite remaining uncertainties, our results suggest that fungal competition may be tightly coupled to soil C storage and sensitive to warming.

We studied how soil fungal diversity, community composition, and decomposition genes responded to 10 years of soil warming, N additions, and their combination. Our results show that soil fungal community composition shifts with warming to a greater extent than N additions, whereas both global changes reduced fungal community evenness by approximately 40%. Our results indicate that global changes can encourage hyperdominance by a handful of mycorrhizal fungal species. Fungal functional potential for the production of manganese peroxidases were reduced while hydrolytic enzyme potentials increased in the heated plots but not in the N addition plots. Warming also reduced negative fungal co-occurrence and promoted species coexistence. Negative co-occurrence was positively correlated with soil organic C content. Through the lens of life history trade-offs, we hypothesize that warming relaxes fungal competition and promotes soil C loss via increased hydrolytic enzyme mediated decomposition. This response by fungi to warming might contribute to a decline in soil C during this phase of soil warming.

All data and code to replicate data manipulation, statistical analyses, and figure generation are available at https://gitlab.ethz.ch/manthony/swan-fungal-responses.

The supplemental files for this article can be found as follows:

Tables S1–S10. Figures S1–S4.

We thank Amber Kittle, Christina Lyons, and Julia Wheeler for laboratory and field assistance. Sequencing was performed by David Miller at the Center for Genomics and Bioinformatics at Indiana University.

This work was funded by a U.S. Department of Defense Strategic Environmental Research and Development Program Grant (NRC2326) to SDF. The Soil Warming × Nitrogen Addition Study at Harvard Forest is maintained with support from the National Science Foundation (NSF) Long Term Ecological Research Program (DEB-1832110) and a Long Term Research in Environmental Biology grant (DEB-1456610) to SDF. MAA was supported by an NSF Graduate Research Fellowship (DGE 1450271) and a Dissertation Year Fellowship from the University of New Hampshire. JAMM acknowledges support from Postdoctoral Development Funds from Oak Ridge National Laboratory. MJS acknowledges support from the Natural Sciences and Engineering Research Council of Canada via a Discovery Grant and a Tier 1 Canada Research Chair in Integrative Molecular Biogeochemistry. Views, opinions, and/or findings contained in this report are those of the authors and should not be construed as an official Department of Defense position or decision unless so designated by other official documentation.

The authors declare no competing interests.

Conceptualization: MAA, SDF.

Formal analysis: MAA, JAMM, MS.

Investigation: MAA, SDF.

Writing—original draft, review, and editing: MAA, SDF, MK, JAMM, MJS.

Supervision: SDF.

How to cite this article: Anthony, MA, Knorr, M, Moore, JAM, Simpson, M, Frey, SD. 2021. Fungal community and functional responses to soil warming are greater than for soil nitrogen enrichment. Elementa: Science of the Anthropocene 9(1). DOI: https://doi.org/10.1525/elementa.2021.000059

Domain Editor-in-Chief: Steven Allison, University of California, Irvine, CA, USA

Associate Editor: Danielle Ignace, Smith College, Northampton, MA, USA

Knowledge Domain: Ecology and Earth Systems

Part of an Elementa Special Feature: The World of Underground Ecology in a Changing Environment