Fungal community response to long-term soil warming with potential implications for soil carbon dynamics

The direction and magnitude of climate warming effects on ecosystem processes such as carbon cycling remain uncertain. Soil fungi are central to these processes due to their roles as decomposers of soil organic matter, as mycorrhizal symbionts, and as determinants of plant diversity. Yet despite their importance to ecosystem functioning, we lack a clear understanding of the long-term response of soil fungal communities to warming. Toward this goal, we characterized soil fungal communities in two replicated soil warming experiments at the Harvard Forest (Petersham, Massachusetts, USA) which had experienced 5°C above ambient soil temperatures for 5 and 20 yr at the time of sampling. We assessed fungal diversity and community composition by sequencing the ITS2 region of rDNA using Illumina technology, along with soil C concentrations and chemistry. Three main findings emerged: (1) long-, but not short-term warming resulted in compositional shifts in the soil fungal community, particularly in the saprotrophic and unknown components of the community; (2) soil C concentrations and the total C stored in the organic horizon declined in response to both short(5 yr) and long-term (20 yr) warming; and (3) following longterm warming, shifts in fungal guild relative abundances were associated with substantial changes in soil organic matter chemistry, particularly the relative abundance of lignin. Taken together, our results suggest that shifts with warming in the relative abundance of fungal functional groups and dominant fungal taxa are related to observed losses in total soil C.


INTRODUCTION
Global warming poses major challenges to the health and functioning of forests worldwide (Peters 1990, Van der Putten 2012, Langley and Hungate 2014, Classen et al. 2015, ecosystems that hold a substantial proportion of carbon (C) in both their vegetation and soil pools (Lal 2004, Bonan 2008, Gill and Finzi 2016. A number of soil warming studies have examined ecosystemlevel effects, with an emphasis on soil respiration and N cycling dynamics (Bradford et al. 2008, Carey et al. 2016, Dacal et al. 2019). Long-term (>10 yr) studies have revealed that chronic soil warming results in a cyclic pattern of soil organic matter (SOM) decomposition and carbon dioxide (CO 2 ) loss to the atmosphere, with phases of substantial soil C loss alternating with phases of no detectable loss (Melillo et al. 2002. Several mechanisms have been proposed to explain this observation including depletion of microbially accessible carbon pools, reductions in microbial biomass, a shift in microbial carbon use efficiency, and changes in microbial community composition (Oechel et al. 2000, Allison and Martiny 2008, Bradford 2013. Of these, the hypothesis that the loss in soil C from long-term warming is driven by changes in microbial community structure is an important one (Schimel and Schaeffer 2012), since warming may selectively impact microbial taxa with specific functional capacities that influence SOM decomposition and, ultimately, rates of soil C cycling.
Soil fungi, integral components of forest ecosystems, mediate key ecosystem processes including C and nutrient cycling (Žifčáková et al. 2016(Žifčáková et al. , Frąc et al. 2018. They form diverse and spatially variable communities as mycorrhizal symbionts, saprotrophs, and pathotrophs (Tedersoo et al. 2014). Fungal community responses to soil warming differ among ecosystems , Treseder et al. 2016); however, a general finding is that the abundance of soil fungi often declines in response to long-term warming (Rinnan et al. 2007, Liang et al. 2015 with shifts in total fungal community composition often observed over the same time frame (Deslippe et al. 2012). It is less clear, however, whether certain functional groups (i.e., ectomycorrhizal vs. saprotrophs) are differentially sensitive to shifts in temperature. For example, the response of ectomycorrhizal fungi to warming has varied from positive (Treseder et al. 2016, Solly et al. 2017 to negative , whereas saprotrophs have been observed to increase in their relative abundance to warming (Treseder et al. 2016, Romero-Olivares et al. 2017. Shifts in abundance of key fungal functional groups (e.g., mycorrhizal fungi vs. saprotrophs; or cellulose vs. lignin-degraders) are expected to significantly impact soil C cycling (Trumbore et al. 1996, Treseder et al. 2016, highlighting the importance of resolving fungal responses to warming.
Here, we build on previous research in a temperate deciduous forest  to assess changes in the diversity and community composition of soil fungi following soil warming. We not only wanted to investigate overall community changes, but also evaluate how different taxa and key functional groups vary in their sensitivities to warming. Since microbial responses to warming have been previously shown to be associated with changes in soil chemistry (Pold et al. 2016, a secondary objective was to assess key soil chemical properties following short-vs. long-term warming as well as examine relationships between the fungal community and soil properties.

Site description, experimental design, and sample collection
This research was performed at two adjacent, replicated soil warming experiments at the Harvard Forest Long-term Ecological Research (LTER) site in Petersham, Massachusetts (MA), USA (42°50 0 N, 72°18 0 W). The experimental plots were established in an area of forest abandoned from agriculture in the late 1800s (Peterjohn et al. 1994). Secondary forest regrowth is dominated by paper and black birch (Betula papyrifera and lenta), red maple (Acer rubrum), black and red oak (Quercus velutina and rubra), and American beech (Fagus grandifolia). Soils are classified as coarse-loamy Inceptisols with no significant variation in pH of the organic soil horizon (Àx = 3.9 AE 0.1 standard error) among experimental plots. Mean annual temperature at Harvard Forest ranges from 19°C to −5°C, with precipitation (annual average~120 cm) distributed evenly throughout the year. One of the longest running soil warming experiments in the world, the Prospect Hill Soil Warming Study, was started in 1991 and consists of six 6 × 6 m control and warmed plots (Melillo et al. 2002). A second experiment, the Soil Warming × Nitrogen Addition Study (SWaN), was established in 2006 adjacent to the Prospect Hill experiment and has six 3 × 3 m control and warmed plots (Contosta et al. 2015). Although there is a nitrogen (N) fertility component at this second experiment, only control and warmed plots were evaluated here. Soils in the warmed plots at both sites are heated continuously to 5°C above ambient temperatures observed in the control plots by the use of resistance heating cables buried 10 cm deep and spaced 20 cm apart .
Soil samples were collected from the organic horizon of four replicate control and warmed plots following five (SWaN) and 20 (Prospect Hill) years of warming for a total of 16 samples (2 sites × 2 warming treatments × 4 replicates). We focused on the organic horizon because previous work at this site has shown that most of the fungal biomass (96%) and microbial activity based on substrate-induced respiration (88%) are located in this horizon , with a significant warming-induced decline in the soil C pool (−31%) being observed for this horizon as well . Samples were collected from experimental plots as intact 20 × 20 cm blocks to the depth of the mineral soil. Samples were sieved (<2 mm), and roots, rocks, and coarse woody debris (>2 mm) were removed. Subsamples were flash-frozen using dry-ice immediately in the field and stored at −80°C until analysis for fungal community characterization. Remaining soils were stored at 4°C until analysis for soil chemical properties.

Fungal community characterization
To characterize fungal diversity and community composition, genomic DNA was isolated on subsamples (0.25 g) using the PowerSoil DNA Isolation Kit (MoBio Laboratories, Carlsbad, California, USA) following manufacturer protocols. Library preparation and sequencing were done at the Department of Energy Joint Genome Institute (Walnut Creek, California, USA). Briefly, polymerase chain reaction (PCR) was performed to amplify the internal transcribed spacer (ITS) 2 region of nuclear rDNA using primers fITS7 (GTGARTCATCGAATCTTTG; Ihrmark et al. 2012) and ITS4 (TCCTCCGCTTATTGATATGC; White et al. 1990). PCRs contained 13 µL of PCRgrade water, 10 µL of Five Prime Hot Master Mix (5 PRIME, Gaithersburg, Maryland, USA), 0.5 µL of each forward and reverse primer (10 µmol/L final concentration), and 1 µL of template DNA. Thermocycling conditions were as follows: one cycle at 94°C for 3 min; then 35 cycles at 94°C for 45 s, 50°C for 60 s, and 72°C for 90 s; and ending with one cycle of 72°C at 10 min. Successful amplification was confirmed using gel electrophoresis (1.5% agarose gel). PCR products were cleaned using the AxyPrep MAG PCR clean-up kit (Coming, Tewksbury, Massachusetts, USA), quantified fluorescently on a Qubit fluorometer (ThermoFisher Scientific, Grand Island, New York, USA), and pooled into equimolar concentrations. Amplicon library sequencing was performed on an Illumina MiSeq v3 (2 × 300 bp). Raw sequences are available at the NCBI database under BioProject number PRJNA523054.

Bioinformatic analysis
Illumina MiSeq data were first demultiplexed and quality-filtered using Trimmomatic v0.36 (Bolger et al. 2014), removing Illumina adapters, sequences with <200 bp, and bases with minimum quality scores <2. All remaining forward and reverse reads were merged with the fastqjoin method using the join_paired_ends.py script in Qiime v1.9 (Caporaso et al. 2010) followed by ITS2 region extraction using ITSx v1.0.11 (Bengtsson-Palme et al. 2013). Sequences were clustered into operational taxonomic units (OTUs) at 97% sequence similarity while simultaneously removing chimeric sequences using the UPARSE-OTU algorithm (Edgar 2013) with the cluster_otus command in USEARCH (v9.2.64; Edgar 2010). Additionally, we excluded all global singletons to reduce artificially inflating richness due to sequencing error. Sequence data were first rarefied to a depth of 257,520 sequences per sample to account for uneven sequence depths across samples. OTUs were taxonomically identified by searching representative sequences of each OTU against the UNITE fungal ITS database (Nilsson et al. 2019) using the BLAST option in the assign_taxonomy.py script in Qiime v1.9. Nonfungal OTUs were excluded from further analyses. Fungal OTUs were initially assigned to three groups (ectomycorrhizal, saprotrophic, and pathotrophic) based on their genus affiliation, trophic mode, and functional guild as previously described (Branco et al. 2013, Tedersoo and Smith 2013, Tedersoo et al. 2014) and using the FUNGuild database (Nguyen et al. 2016). OTUs were placed into groupings only if assignments were deemed as highly probable or probable based on default parameters in the FunGuild v www.esajournals.org database (http://stbates.org/guild/app.php). We further divided ectomycorrhizal fungal OTUs into ectomycorrhizal fungal exploration types (i.e., contact to short distance or medium to long distance) according to (Agerer 2001), while saprotrophic fungal OTUs were also further divided into yeasts, filamentous saprotrophs, and white rot fungi according to the FUNGuild database (Nguyen et al. 2016). Here, white rot fungi were included as a separate group because of our interest in lignin-degraders. Fungal OTUs were designated as "taxonomy insufficient" if they were assigned a genus affiliation but their functional grouping was uncertain. Fungal OTUs that either had no known identification in available databases or were identified only to the kingdom level were designated as "unknown taxonomy."

Soil chemical properties
Soil bulk density, determined for each sample, was used to calculate volumetric soil C stocks (g C/m 2 ). Soil samples were analyzed for total organic C and N on air-dried finely ground samples by dry combustion on a Perkin Elmer 2400 Series II CHN elemental analyzer (Waltham, Massachusetts, USA). Organic matter chemistry was assessed by determining the relative abundances of chemical classes using pyrolysis-GC/ MS (Wickings et al. 2011). In brief, samples were pyrolyzed at 600°C, products transferred to a GC where compounds were separated and transferred to an ion trap mass spectrometer, ionized, detected via electron multiplier, and identified using a compound library from the National Institute of Standards and Technology (NIST; Grandy et al. 2007Grandy et al. , 2009).

Statistical analyses
All statistical analyses were carried out using R v3.3.2 (R Core Team 2019). Due to the close physical proximity of the experimental plots for the two warming experiments (Prospect Hill and SWaN), all control plots were grouped together for downstream analyses as no significant differences were detected across control plots for any measured soil property, including fungal community structure (Appendix S1: Table S1). We ran linear models to examine whether there were differences in measured soil chemical properties in response to short-(5 yr) or long-term (20 yr) warming. Generalized linear models with a Poisson distribution were run to examine how fungal OTU richness responded to warming, while linear models were run to examine whether there were differences in the relative abundance of soil fungal taxa or functional groups to warming. Post hoc pairwise comparisons were performed using the glht function in the package multcomp (Hothorn et al. 2008). All model assumptions were checked by visual inspection of residual diagnostic plots (Zuur et al. 2009).
We used multivariate generalized linear models (M-GLM) to test whether fungal OTU relative abundance distributions were affected by warming with the manyglm function in the mvabund package (Wang et al. 2012), using a negative binomial model, likelihood ratio statistic, and P values assigned following 9999 iterations via probability integral transform-trap sampling to account for correlation in testing using the anova function. We used this approach instead of a distance-based metric such as perMANOVA, as mvabund directly accounts for the strong mean-variance relationships present in most multivariate community relative abundance analyses (Warton et al. 2015). To visualize differences in soil fungal communities across warming treatments, we used non-metric multidimensional scaling (NMDS) ordination based on Bray-Curtis dissimilarity using the metaMDS function in the vegan package (Oksanen et al. 2014). To identify fungal OTUs common or representative of either control or warmed plots, we used the multinomial species classification method in the vegan package (Oksanen et al. 2014) with supermajority rule and all other parameters set as default (i.e., CLAM test) to statistically classify fungal OTUs into the following categories: fungi shared across control and warmed plots, fungi preferentially associated with either control or warmed plots, and fungi that were too rare (based on low relative abundance) to be assigned association specificity (Chazdon et al. 2011).
Partial least-square regression (Wold et al. 2001) was used to identify covariates (soil chemistry) most strongly correlated with fungal guild richness and relative abundance using the plsr() function in the pls package (Mevik and Wehrens 2007). Models were refined to the most important predictor variable(s) and analyzed independently using linear regression. All model v www.esajournals.org assumptions were checked by visual inspection of residual diagnostic plots (Zuur et al. 2009).

Fungal diversity and community composition
In total, 6,988,799 sequences were obtained across all samples. After quality filtering, 6,165,723 sequences (88%) representing 2552 fungal OTUs were assessed for taxonomic affiliation. At the phylum level, there were 1094 OTUs (43%) assigned to Basidiomycota, 498 OTUs (19%) assigned to Ascomycota, and 107 OTUs (4%) assigned to basal clades, while Glomeromycota made up <1% of identified sequences (Appendix S1: Table S2). Since ITS primers are biased against glomeromycetes (arbuscular mycorrhizal fungi; Stockinger et al. 2010), this group was not included in subsequent analyses. At the functional group level, there were 732 OTUs assigned to ectomycorrhizal fungi (28%), 339 OTUs assigned to saprotrophic fungi (13%), and 36 OTUs assigned to pathotrophic fungi (1% ;  Table 1). In addition, there were 843 fungal OTUs (33%) assigned as unknown taxonomy, with most not being assigned to a phylum. A final group of 593 (23%) OTUs were assigned a taxonomic affiliation; however, there was insufficient ecological information with which to reliably assign them to a functional group.
Overall, long-(20 yr) but not short-term (5 yr) warming resulted in significantly different soil fungal communities compared to the control, and long-term plots were also significantly different from short-term warming plots ( Fig. 1; Appendix S1: Table S3). Fungal OTU richness tended to increase with warming compared to control plots; however, this increase was only significant for saprotrophic fungi (Table 2). Specifically, the richness of filamentous saprotrophs increased from 61 OTUs in the control treatment to 80 OTUs following 20 yr of warming. In contrast, saprotrophs declined in relative abundance from 33.7% in control plots to 11.6% in plots warmed for 20 yr. This decline was largely due to a significant negative response to warming of filamentous saprotrophs, which made up the majority of saprotrophic fungi (~74%). As the relative abundance of saprotrophic taxa declined with long-term warming, that of unknown fungi significantly increased, going from 28.7% in the control treatment to 41.5% in warmed soils. Ectomycorrhizal and pathotrophic fungal communities did not differ compositionally or in their relative abundances in control vs. warmed plots (Fig. 1, Table 2). Although fungal biomass was not measured as part of this study, previous research at our site has shown that fungal biomass is not significantly affected by short-term warming (<5 yr; Contosta et al. 2015), while fungal biomass tends to decline in response to warming after 12-20 yr .
Since most significant differences were observed in response to long-term warming, we further investigated whether there were representative and/or unique taxa in soils warmed for 20 yr. Overall, about a third of fungal OTUs were present in both control and warmed (20 yr) soils (Appendix S1: Fig. S1a), whereas about half of all OTUs were found primarily in either one or the Notes: Values represent the total of all OTUs across experimental treatments. Values in parentheses are the percentages represented by each category based on total assigned OTUs. OTUs that could not be identified to the family level were assigned a function of unassigned. General categories are bolded for clarity.
† Represents fungal OTUs that were assigned a taxonomic affiliation but for which their functional guild was uncertain or unknown.
‡ Includes OTUs that were identified to the kingdom level (i.e., fungi), but which have no known taxonomic identification in available databases. other treatment, with an equal proportion observed to be unique in control (340 OTUs) compared to warmed (345 OTUs) soils. In terms of specific functional groups, a third of ectomycorrhizal fungal OTUs were present in both treatments (Appendix S1: Fig. S1b), while saprotrophic fungi and unknown fungi were found primarily in warmed (saprotrophic, 31%, 65 OTUs; unknown, 29%, 140 OTUs) compared to control soils (saprotrophic, 15%, 31 OTUs;  Note: Values represent means and standard errors, with different lowercase superscript letters in bold indicating significant difference (P < 0.05).
† Represents fungal OTUs that were assigned a taxonomic affiliation but for which their functional guild was uncertain or unknown.

Soil chemical properties
Warming resulted in a significant decline in soil C concentrations (F = 3.72, P = 0.04) and in the total C stored in the organic horizon (F = 5.46, P = 0.01) over short-(5 yr) and longterm (20 yr) warming (Table 4). Soil organic matter chemistry was also altered, but only following long-term (20 yr) warming. More specifically, two decades of warming induced a significant reduction (−37%) in the relative abundance of lignin. In contrast, the relative abundance of polysaccharides and unknown origin chemical classes were significantly higher over the same warming period. ‡ Represents fungal OTUs that were assigned a taxonomic affiliation but for which their functional guild was uncertain or unknown.

Relationships between fungi and soil chemistry
Overall, there were no significant soil chemical predictors of fungal guild richness. However, the relative abundance of lignin was found to be significantly correlated to fungal guild relative abundance (Fig. 2). The relative abundance of saprotrophic fungi as a whole and the most abundant saprotroph, Megacollybia, were positively correlated with lignin (saprotrophic fungi, F = 15.95, P = 0.001; Megacollybia, F = 13.92, P = 0.002). The relative abundances of ectomycorrhizal and unknown fungi were negatively correlated with lignin (ectomycorrhizal: F = 9.26, P = 0.009; unknown: F = 3.49, P = 0.05). In contrast, there were no significant correlations with pathotrophic relative abundance.

DISCUSSION
We assessed soil fungal community diversity and composition, soil C stocks, and key soil chemical properties following short-and longterm warming. Three main findings emerged: (1) long-, but not short-term warming resulted in compositional shifts in the soil fungal community; (2) soil C concentrations in the organic horizon significantly declined in response to short-(5 yr) and long-term (20 yr) warming; and (3) following long-term warming, shifts in fungal guild relative abundances were associated with substantial changes in soil organic matter chemistry, particularly a decline in the relative abundance of lignin.
We found that the composition of the overall soil fungal community was significantly different in warmed vs. control plots, but only following long-term (20 yr) warming. This community shift appears to be attributable to certain fungal functional guilds being differentially sensitive to warming. For example, saprotrophic fungi have been shown to generally increase in relative abundance with warming (Treseder et al. 2016, Romero-Olivares et al. 2017, while the responses of ectomycorrhizal fungi are varied, being dependent on their hosts' growth and photosynthetic rates (Clemmensen et al. 2006, Fernandez et al. 2017. Our contrasting findings of decreased relative abundance of saprotrophic fungi with soil warming are potentially due to (1) competition with ectomycorrhizal fungi (Averill et al. 2014) as supported by an increased ratio of ectomycorrhizal to saprotrophic fungi in long-term warmed plots (2.33 AE 0.2) compared to control plots (0.70 AE 0.5), and/or (2) lower C concentrations and lignin relative abundance. There were shifts in the most abundant soil fungal taxa found in response to 20 yr of soil warming. For example, Megacollybia, a common saprotroph in hardwood and mixed forests of northeastern North America, is a highly competitive basidiomycete shown to dominate decomposition (Boddy and Watkinson 1995). It is a white rot fungus with lignin-degrading capability (Worrall et al. 1997). In our study, the relative abundance of Megacollybia was significantly correlated with lignin, declining threefold in warmed vs. control plots (Fig. 2D). We hypothesize that as the relative abundance of lignin declines in response to warming, so does that of Megacollybia. While the relative abundance of Megacollybia declined, that of Oidiodendron, a dark septate endophyte with a saprotrophic lifestyle (Nilsson 1973) increased in response to warming. Oidiodendron is capable of producing a range of enzymes necessary for cellulose and lignin decomposition (Haselwandter et al. 1990, Rice andCurrah 2005) and has been observed previously to respond positively to warming (Solly et al. 2017).
Surprisingly, we found ectomycorrhizal fungi to be invariant to soil warming. Though global warming can alter the composition of ectomycorrhizal fungi by pushing them and/or their host plants outside their natural range of physiological tolerances (Kipfer et al. 2010, Fernandez et al. 2017, several other studies examining ectomycorrhizal response to increased temperatures have found small and/or absent effects (Parrent et al. 2006, Tu et al. 2015, Fernandez et al. 2017, Mucha et al. 2018. For example, in the boreal zone, ectomycorrhizal colonized tree species exposed to warming have shown a reduction in growth and photosynthetic rate , Fernandez et al. 2017, while warming of tree species in temperate zones has not resulted in similar reductions Oleksyn 2008, Wheeler et al. 2017). This suggests that in regions where host plant performance declines with increased temperature, trees may allocate less C to ectomycorrhizal fungi, resulting in declines in ectomycorrhizal fungal richness and relative abundance.
Fungal taxa of unknown identity showed a strong positive response to long-term warming, significantly increasing in their richness and relative abundance compared to the control treatment. Unknown fungi also showed a significant negative correlation with soil chemistry, where lignin relative abundance correlated with lower relative abundance of unknown taxa. In addition, a number of taxa with unknown identity were unique indicator species of either the control treatment or plots warmed for 20 yr. Identifying and understanding the specific v www.esajournals.org contributions of the dominant unknown fungal taxa described here are of high priority as they show a strong response to long-term warming and may play a critical role in soil C dynamics. The proportion of unidentified taxa (33% of OTUs) is consistent with other fungal community analyses employing high-throughput sequencing (Lücking and Hawksworth 2018 and references therein). The poor characterization of fungi is due to a lack of reference sequences and mis-annotated or misidentified reference sequences (Nilsson et al. 2011, Lücking and Hawksworth 2018. Initiatives such as the UNITE database for the molecular identification of fungi (Koljalg et al. 2013) and alternative classification and identification methods (Hibbett et al. 2016, Irannia and Chen 2016, Abarenkov et al. 2018, Hongsanan et al. 2018) have been proposed to mitigate these issues.
Though this current study focused on fungi, the importance of bacteria in our experimental system at Harvard Forest has been separately studied. Interestingly, bacteria showed a similar response to warming, where only the long-term warming (20 yr) altered bacterial community composition . The fungalto-bacterial ratio of the whole microbial community was also reduced in the warming plots caused by a reduction in fungi, not bacteria , and warming also led to an enrichment of putative lignin-using bacterial taxa . This suggests that the decline in relative abundance of lignin may be facilitated by microbial community shifts toward bacterial taxa with the capacity to enzymatically access and use lignin, leading to increased competition between bacteria-and fungal lignin-degrading taxa (e.g., Megacollybia).
In summary, long-term warming induced changes in fungal diversity and community composition, as well as the relative abundance of fungal groups with different trophic strategies. Our results also highlight that long-term warming has rather large net effects on soil C storage and soil organic matter chemistry, particularly the relative abundance of lignin. Potential reductions in soil C storage and organic matter quality and quantity (Kammer et al. 2009, Harte et al. 2015 could, in turn, effect future ecosystem C fluxes and storage. Taken together, we propose that shifts with warming in the relative abundance of key fungal functional groups and dominant fungal taxa, along with lower total fungal biomass , are related to losses in soil C concentrations and the total C stored in the organic horizon.

ACKNOWLEDGMENTS
Thanks to Brian Godbois and Chris Cook at the University of New Hampshire for soil sample collection. Sequencing was done at the Joint Genome Institute as part of a Community Sequencing Program (under contract no. DE-AC02-05CH11231 CSP-1058) award. Thanks to Toni Westbrook at UNH for assistance with bioinformatics. The soil warming plots at Harvard Forest are maintained through support from the NSF Long-term Ecological Research (DEB 1237491) and Long-term Research in Environmental Biology (DEB 1456528) Programs.