Immediate and legacy effects of snow exclusion on soil fungal diversity and community composition

Soil fungi play crucial roles in ecosystem functions. However, how snow cover change associated with winter warming affects soil fungal communities remains unclear in the Tibetan forest. We conducted a snow manipulation experiment to explore immediate and legacy effects of snow exclusion on soil fungal community diversity and composition in a spruce forest on the eastern Tibetan Plateau. Soil fungal communities were performed by the high throughput sequencing of gene-fragments. Ascomycota and Basidiomycota were the two dominant fungal phyla and Archaeorhizomyces, Aspergillus and Amanita were the three most common genera across seasons and snow manipulations. Snow exclusion did not affect the diversity and structure of soil fungal community in both snow-covered and snow-free seasons. However, the relative abundance of some fungal communities was different among seasons. Soil fungal groups were correlated with environmental factors (i.e., temperature and moisture) and soil biochemical variables (i.e., ammonium and enzyme). These results suggest that the season-driven variations had stronger impacts on soil fungal community than short-term snow cover change. Such findings may have important implications for soil microbial processes in Tibetan forests experiencing significant decreases in snowfall.


Background
Seasonal snow cover is one of the most important factors that drive biogeochemical cycling in cold regions (Jusselme et al. 2016;Liu et al. 2018). Winter warming is predicted to reduce the stability and thickness of snow cover in snowy regions (Gobiet et al. 2014). The decline or absence of insulation of the snow cover would alter soil environmental conditions (e.g., temperature, moisture, frost intensity and duration) during winter and the growing season (Wipf and Rixen 2010;Kreyling et al. 2012;Aanderud et al. 2013;Li et al. 2017;Song et al. 2017). Soil environmental conditions are very important to mediate soil microbial community in cold ecosystems (Zinger et al. 2009;Voříšková et al. 2014;Morgado et al. 2016). Hence, reduced snow cover associated with winter climate change can result in more intensive soil freezing, which may in turn have profound influences on soil biological processes, especially soil fungal communities.
Soil fungi are the key decomposer of soil organic matter in cold soils (Baldrian 2017). The diversity and composition of fungal communities play key roles in soil carbon and nutrient cycling (Cheng et al. 2017;Asemaninejad et al. 2018). Therefore, it is crucial to understand how fungal community responds to snow cover change in alpine soils. Some studies have shown that snow depth change resulted in an immediate effect on soil fungal community composition in snow-covered winter (Olofsson et al. 2011;Barbeito et al. 2013;Voříšková et al. 2014;Santalahti et al. 2016). However, winter snowpack can affect the soil moisture and nutrients in the subsequent growing season (Wipf and Rixen 2010). Therefore, snow cover change could also lead to carryover effects on soil fungal community composition in snow-free growing season (Buckeridge et al. 2013;Wubs et al. 2018;Sorensen et al. 2020). To our knowledge, soil fungal responses to altered snow cover have scarcely been investigated both in snow-covered winter and snow-free growing season. Obviously, exploring the immediate and legacy effects of snow cover change on fungal community is very essential to understand microbe-associated ecological processes in cold soils.
As the earth's 'Third Pole', the Tibetan Plateau is anticipated to become warmer in the coming decades (Chen et al. 2013). In this region, winter snowfall has been decreasing over last decades (Wang et al. 2016;Deng et al. 2017). Winter soil temperature is close to the physical melting point and is sensitive to snow cover change Li et al. 2017). Our prior study has found that snow exclusion reduced soil respiration and enzyme activities in wintertime, but did not result in any cross-seasonal effects in subsequent snowfree season (Yang et al. 2019). Fungal communities are considered as key agents controlling soil C cycling in cold ecosystems (Clemmensen et al. 2013;Zhang et al. 2013;Solly et al. 2017). In this study, we aimed to evaluate the immediate and legacy effects of snow exclusion on the diversity and composition of soil fungal communities. Specifically, we hypothesized that, snow exclusion (1) will decline soil fungal communities, and thereby change fungal community composition; (2) will result in both immediate and carry-over impacts on soil fungal communities; (3) will affect fungal communities via changes in environmental and/or soil biochemical factors.

Site description
The study was conducted in a dragon spruce (Picea asperata) stand at the Long-term Research Station of Alpine Forest Ecosystems of Sichuan Agricultural University, which is located at the eastern Tibetan Plateau of China (31°15′ N, 102°53′ E; 3021 m a.s.l.). The mean annual temperature is 3.0°C, with maximum and minimum temperatures of 23.0°C (July) and − 18.0°C (January), respectively. Annual precipitation is about 850 mm. In general, snow begins to accumulate in late November and melts in late March of the following year. The understory is dominated by Salix paraplesia, Rhododendron lapponicum, Cacalia sp., Carex sp., and Cyperus sp. (Li et al. 2017). The soil is classified as Cambic Umbrisols (IUSS Working Group WRB 2007).

Experimental design
Winter snowfall was excluded using shelters. This technique can effectively reduce snow cover and minimize unwanted environmental conditions (Li et al. 2017). In early November 2015, six wooden roofs (2 m height, 3 m × 3 m ground area) were established in the Picea asperata forest to prevent snow accumulation on the ground. One control plot that allows snow input was set up in the vicinity of each wooden roof (3 m × 3 m ground area). In late winter, the accumulated snow on the roof was added to the forest floor in order to ensure the similar water balance between the snow-free and control plots. The snow manipulation began in late November 2015 and ended in early April 2016 when the seasonal snow in the control plots was melted (Li et al. 2017).

Microclimate and soil biochemical analyses
Soil temperature (5 cm depth) and air temperature (2 m height) were measured every 1 h by the Thermochron iButton DS1923-F5 Recorders (Maxim Dallas Semiconductor Corp., USA). The minimum daily mean soil temperatures of the control and snow exclusion plots were − 0.5°C and − 2.2°C, respectively (Fig. S1). The snow exclusion increased winter soil frost and lowered the average soil temperature (Fig. S1a). Seasonal snow began to accumulate in late November 2015 and melted in early April 2016 (Fig. S1). Snow depth in the control plots was measured approximately every 2 weeks. There was no obvious difference in soil moisture between control and snow exclusion plots (Fig. S1b).
Soil pH was determined on field moist soil in a 1:2.5 (M/V) soil suspension using a pH meter (PHS-25CW, BANTE Instruments Limited, Shanghai, China). Soil ammonium ( NH þ 4 -N) and nitrate ( NO − 3 -N) were extracted with 2 mol•L − 1 KCl, and then NH þ 4 and NO − 3 in the extracts was determined using colorimetry (Xu et al. 2010). Soil microbial biomass carbon (MBC) was measured by the fumigation-extraction method (Vance et al. 1987). The released C was converted to MBC using kec 0.45 (Vance et al. 1987).

Soil sampling
In 2016, three snow exclusion plots and their corresponding controls were randomly selected. Soil samples were collected from the topsoil (0-15 cm) in mid-February (deep snow period, DSP), early April (early thawing period, ETP), and mid-August (middle of the growing season, MGS), respectively. On each sampling date, three soil cores (5 cm in diameter, 0-15 cm deep) were randomly taken at each plot and were mixed into one composite sample per plot. The composite sample was passed through a 2-mm sieve, and any visible living plant material was removed from the sieved soil. Subsamples of the sieved soils were stored in the refrigerator at − 70°C and 4°C, respectively, for DNA and routine chemical analyses.

DNA extraction, PCR amplification and Illumina MiSeq sequencing
Microbial DNA was extracted from 18 soil samples by using the E.Z.N.A.® Soil DNA Kit (Omega Bio Inc. Norcross, GA, USA) according to manufacturer's protocols. The final DNA concentration and purification were determined by Nanodrop® ND-1000 UV-Vis spectrophotometer (Nano-Drop Technologies, Wilmington, DE, USA), and DNA quality was checked by 1% agarose gel electrophoresis.
The fungal 18S rDNA gene was amplified with the primers SSU0817F/SSU1196R (Borneman and Hartin 2000) by a thermocycler PCR system (GeneAmp 9700, ABI, USA). PCR reactions were performed in triplicate using a 20 μL mixture containing 4 μL of 5 × FastPfu Buffer, 2 μL of 2.5 mmol•L − 1 dNTPs, 0.8 μL of each primer (5 μmol•L − 1 ), 0.4 μL of FastPfu Polymerase and 10 ng of template DNA. The PCR reactions were conducted using the following program: denaturation at 95°C for 3 min, annealing at 55°C for 30 s, and elongation at 72°C for 45 s, and a final extension at 72°C for 10 min. The step from the denaturation to extension was run for 35 cycles.
The resulted PCR products were extracted from a 2% agarose gel and further purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) and quantified using QuantiFluor™-ST (Promega, USA) according to the manufacturer's protocol. The purified amplicons were pooled in equimolar and paired-end sequenced (2 × 300) on an Illumina MiSeq platform (Illumina, San Diego, USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw reads were deposited in the NCBI Sequence Read Archive (SRA) database with accession number PRJNA562976.

Processing of sequencing data
Raw FASTQ files were demultiplexed and qualityfiltered by Trimmomatic and merged by FLASH with the following criteria: (i) the reads were truncated at any site receiving an average quality score < 20 over a 50-bp sliding window; (ii) Primers were exactly matched allowing two-nucleotide mismatching, and reads containing ambiguous bases were removed; (iii) Sequences whose overlap was longer than 10 bp were merged according to their overlap sequence.
All sequences acquired using the Illumina-MiSeq was saved in the raw fastq files. Initial processing of the raw dataset included screening to remove short and lowquality reads; only high-quality sequences were retained. Operational taxonomic units (OTUs) were clustered with 97% similarity cutoff using UPARSE (version 7.1, http:// drive5.com/uparse/), and chimeric sequences were identified and removed using UCHIME (Edgar et al. 2011). The taxonomy of each 18S rDNA gene sequence was analyzed by the RDP Classifier algorithm (http://rdp.cme.msu.edu/) against the Silva 128/18S_eukaryota database using a confidence threshold of 70%. Rarefaction curves and alpha diversity calculations were based on OTUs with > 97% identity. Rarefaction analysis and alpha-diversity indices (abundance-based Sobs, Chao1, Shannon and Simpson) were revealed by Mothur (Schloss et al. 2009).

Statistical analysis
Alpha diversity metrics, including the Shannon-Wiener index, Simpson's diversity index, richness (Sobs) and community coverage, were calculated to determine the "diversity" and "richness" functions of the fungal community. Repeated measures ANOVAs were performed to test the effects of treatment, sampling date (deep snow period, early thawing period, middle of the growing season), and their interactions on the fungal indices. For specific sampling dates, Student t-tests were used to compare the effect of snow exclusion. The statistical tests were considered significant at the P < 0.05 level. All statistical analyses were performed using SPSS 20.0 (IBM Corporation, Armonk, NY, USA).
One-way ANOVA with Tukey-kramer post hoc tests were performed to test the effects of sampling date on the abundance of fungal community at the same treatment. For individual sampling dates, Wilcoxon ranksum test was used to compare the effect of snow exclusion on the abundance of fungal community structure. Statistical testing among variation in fungal community composition was carried out using the analysis of similarity (ANOSIM). Differences in soil fungal phyla were represented on a two-dimensional ordination plot following Non-metric multidimensional scaling (NMDS) analysis based on Bray-Curtis distance. Spearman correlation heatmap analysis was performed to examine the relationships between the relative abundance of fungal taxa and environmental factors or biochemical properties. All statistical analyses were performed using the VEGAN package of the R software (Oksanen et al. 2013;R Development Core Team 2015).

Pyrosequence data description and species alpha diversity
Pyrosequencing of 18 samples had a total of 604,373 raw sequences. A total of 215,046 sequences were retained after quality control. The sequence data were classified into 167 OTUs at 97% similarity. After normalization, each library had 11,947 reads and the sequence was clustered into 70-137 OTUs (Table S1). All rarefaction curves tended to approach the saturation plateau, indicating that the data volume of sequenced reads was reasonable (Fig. 1a). In the deep snow period, the Shannon index was higher compared to the early thawing period and the middle of the growing season (Fig. 1b). However, the Shannon and Simpson indices were unaffected by snow exclusion (Fig. 1b).

Variation in fungal community composition
Across all sites, fungal communities were consistently dominated by Ascomycota (38.3%-55.6%) and Basidiomycota (7.92%-54.50%) (Fig. 2a). The relative abundance of the Chytridiomycota (0.81%-18.1%), Glomeromycota (0.33%-1.18%), LKM15 (0.05%-0.84%) and Cryptomycota (0.01%-0.19%) was very low in all samples. Based on Bray-Curtis, ANOSIM and NMDS were used to compare the similarity of the fungal phyla between treatments or seasons (Fig. 3). There were significant differences in soil fungal communities among seasons. Soil fungal community in the middle of the growing season was obviously different from that in the deep snow period and early thawing period ( Fig. 3; ANOSIM, r = 0.527, P = 0.001). In the control plots, Basidiomycota was much higher in the middle of the growing season than in the deep snow period (P < 0.05) (Fig. S2). In the treatment plots, LKM15 in the deep snow period was much higher than that in the middle of the growing season (P < 0.05) (Fig. S2). However, no significant differences were found between snow regimes (Figs. 3 and S3; ANOSIM, r = − 0.031, P = 0.539).
At the genus level, 17 distinct groups were observed in all samples which belonged to 12 classes (Fig. 2b, Table 1). Archaeorhizomyces (Archaeorhizomycetes), Aspergillus (Eurotiomycetes) and Amanita (Agaricomycetes) were the first three abundant genera (Fig. 2b, Table 1). Snow exclusion did not show any significant effects on the abundance of each genus, while sampling period had significant effect (Figs. S4 and S5).

Correlation between soil biochemical properties and fungal communities
Spearman correlation heatmap analysis was performed to examine the relationships between the phylum and genus level of fungal communities and soil environmental and biochemical variables (Fig. 4). For the known phyla, soil temperature, NH þ 4 -N and enzyme activities (BG, AP and NAG) were important influence factors (Fig. 4a). For the three common genera, a negative correlation was observed between Aspergillus and NO − 3 -N. Moreover, Amanita showed significant correlations with pH, NH þ 4 -N, BG and AP activities (Fig. 4b).

Discussion
The depth and duration of winter snow have been decreasing on the eastern Tibetan Plateau (Li et al. 2017).
Recent studies have reported that snow cover change caused significant effects on soil fungal communities in arctic cold ecosystems Semenova et al. 2016). We utilized a snow-exclusion experiment to test the immediate and cross-seasonal effect of snow exclusion on the diversity and composition of soil fungal communities in a subslpine spruce forest on the eastern Tibetan Plateau. Inconsistent with our hypotheses, snow exclusion did not affect the diversity and composition of soil fungal communities in both snow-covered winter and snow-free growing season, but the fungal diversity and composition greatly varied across seasons, indicating  Fig. 1 Rarefaction curve of the OTU number at 97% similarity cutoff (a) and α diversity index (b) of soil fungal community in control (C) and snow exclusion (SE) plots in the deep snow period (DSP), early thawing period (ETP), and in the middle of the growing season (MGS). * indicated a significant difference between control and snow exclusion on a specific same sampling date. Data shown are mean ± s.e that soil fungal communities were not sensitive to shortterm mild increased frost associated with snow cover absence. Both phylum Ascomycota and Basidiomycota were the two dominant phyla in the study site. Similar findings are observed in other cold ecosystems, such as glacier ecosystems, alpine meadows and boreal forests Antony et al. 2016;Gao and Yang 2016;Chen et al. 2017aChen et al. , 2017bHan et al. 2017;Männistö et al. 2018). Hypocreales, Archaeorhizomycetales and Eurotiales were the dominant orders in the phylum Ascomycota. In addition, Agaricales was the dominant order in the phylum Basidiomycota. As stated above, there might be similar common soil fungal communities among different snowy environments, indicating these fungal communities have an extensive adaptability.   A growing number of studies have revealed that snow cover changes leaded to significant impacts on soil fungal community composition and diversity in cold ecosystems (Mundra et al. 2016;Semenova et al. 2016;Männistö et al. 2018). For example, deeper snow cover decreased saprotrophic fungi but increased ECM fungi richness (Mundra et al. 2016). In addition, increased snow cover significantly altered the composition of soil fungal communities in arctic tundra Semenova et al. 2016). In contrast, the lack of snow cover resulted in slight effects on soil fungal community structure and activity in a boreal coniferous forest (Männistö et al. 2018). Our results are consistent with the observations found in alpine tundra (Zinger et al. 2009) and in temperate and boreal forest ecosystems (Gao et al. 2018), which indicated snow cover change did not affect soil fungal diversity and community composition. This is probably because soil fungal communities in Fig. 3 Non-metric multidimensional scaling (NMDS) plot showing variation in the composition (Bray-Curis distance) of soil fungal communities between snow regimes and sampling periods on phylum level. Different symbols represent different sampling date. The inverted triangle, circle and diamond represents the deep snow period (DSP), early thawing period (ETP) and in the middle of the growing season (MGS), respectively. The hollow represents Control (C) and solid represents snow exclusion (SE) plots Fig. 4 Correlation heatmap of soil biochemical properties and read numbers at the phylum (a) and genus (b) level for fungal. NH þ 4 -N: ammonium nitrogen; NO − 3 -N: nitrate nitrogen; MBC: microbial biomass carbon; BG: 1, 4-β-glucosidase; AP: acid phosphatase; NAG: β-N-acetylglucosaminidase. The color intensity in each panel indicates the relative correlation between soil property and read numbers of each group. * P ≤ 0.05, ** P ≤ 0.01, *** P ≤ 0.001 cold regions have developed physiological resistance to extreme conditions, such as low temperature and freezethaw cycles (Stres et al. 2010;Haei et al. 2011). Additionally, fungi can still gain energy through dissolved organic carbon and nitrogen in frozen soils, implying that soil fungal communities had some unique adaptive strategies to survive extreme cold conditions (Fitzhugh et al. 2001;Matzner and Borken 2008;Sorensen et al. 2016). As a result, a short-term mild soil freezing is not powerful enough to alter the diversity and composition of fungal communities in alpine forest soils. On the other hand, snow exclusion did not affect soil fungal communities in the middle of the growing season, implying that the carryover effect of snow exclusion is negligible. This is mainly because no significant differences were detected in both environmental factors (e.g., temperature and moisture) and soil variables (e.g., N pools, enzyme activities) in this period (Yang et al. 2019).
Although no significant treatment effects were detected, fungal community diversity and composition significantly varied across seasons. This observation is consistent with the results found in boreal forest soils (Voříšková et al. 2014). The Shannon-Wiener index in the deep snow period is higher than those in the early thawing period and in the middle of the growing season, which may be attributed to the low plant cover and diversity in winter (Shi et al. 2014). A study conducted in a Scots pine forest showed that saprotrophs are dominant in winter, but ECM fungi grow rapidly in the growing season (Santalahti et al. 2016). In this study, the Basidiomycota group and genera Amanita (Agaricomycetes) showed a distinct transition from frozen winter to the subsequent growing season. This is mainly because the carbon sources for Basidiomycota are mostly derived from exogenous materials, such as plant litter and wood (Kellner et al. 2010). Therefore, relatively large vegetation coverage and abundance may provide more substrates during the growing season. In addition, relatively rich root exudates can also favor soil fungal community, especially in colder ecosystems (Shahzad et al. 2015;Delgado-Baquerizo et al. 2019). In general, Amanita has a mycorrhizal relationship with vascular plants (Yang 2000). Therefore, plant growth in the growing season may stimulate Amanita growth. At the same time, the warmer temperature in the growing season may also provide suitable conditions for Amanita reproduction and growth (Asemaninejad et al. 2017).
Ascomycota did not vary among seasons in this case, indicating that these fungi having strong resistance and adaptability to environmental stress. Some studies have demonstrated that Ascomycota can adapt to harsh habitats to obtain survival advantages (Chen et al. 2017a;Asemaninejad et al. 2018). However, the genus members of Ascomycota, Archaeorhizomyces and Aspergillus shifted across seasons. The relative abundance of Archaeorhizomyces was found to be highest in the growing season. Studies showed that Archaeorhizomyces, as a group of saprophytic fungi, are dominant in summer (Schadt et al. 2003, Fig. 3) and omnipresent in roots and rhizosphere soil (Rosling et al. 2011). Root-derived compounds are the principal carbon source for these fungi (Rosling et al. 2011). Thus, sufficient high-quality substrates associated with root activity are favorable to Archaeorhizomyces in the growing season (Lindahl et al. 2007). Conversely, the relative abundance of Aspergillus was higher in the deep snow period as compared to the other two periods, which is consistent with the observations in a Pal forest soil (Rane and Gandhe 2006).
Environmental and biochemical variables are two important drivers of microbial community composition (Fitzhugh et al. 2001;Matzner and Borken 2008;Chen et al. 2017b;Asemaninejad et al. 2018). Soil temperature is a key factor affecting fungal communities (Asemaninejad et al. 2017). The abundance of Basidiomycota is relatively high in the warm growing season (Kirk et al. 2001;Asemaninejad et al. 2017). In addition, soil moisture also regulates soil fungal community in forest ecosystems (Bainard et al. 2014;Kim et al. 2016). The significant differences in soil moisture between winter and growing season may, to some extent, account for the separation of fungi group, such as Basidiomycota and Agaricomycetes. Soil enzyme production is mainly derived from soil microbes, especially soil fungi in forest ecosystems (Schneider et al. 2012). Thus, soil enzymes can partially reflect the composition and structure of fungal communities (Frey et al. 2004;Kivlin and Treseder 2014). In our study, soil enzyme activities (e.g., BG, AP, and NAG) showed significant correlations with some specific fungal communities. Similar observations were also found in the California forests ecosystem (Kivlin and Treseder 2014). Soil nutrients, especially nitrogen availability, could affect soil microbial community in cold biomes (Bardgett et al. 2002). Many studies have demonstrated that nitrogen additions alter soil fungal community composition (Pardo et al. 2011;Morrison et al. 2016;Corrales et al. 2017). In this case, seasonal variations in N pools may partly explain the seasonal dynamic of fungal communities. Previous studies have demonstrated that soil fungal communities have wider pH ranges for optimal growth (Rousk et al. 2010). Similarly, soil pH was not correlated with fungal communities at the phylum level, suggesting that soil pH is less important to mediate soil fungal communities in subalpine coniferous forests.

Conclusions
This study examined the immediate and legacy effects of snow cover change on the diversity and composition of fungal communities in a subalpine spruce forest on the Tibetan Plateau of China. Our findings suggested that snow exclusion did not affect soil fungal communities in both snow-covered winter and snow-free growing season, indicating that fungal communities were insensitive to short-term snow cover change in Tibetan forest soils. Besides, soil fungal communities varied across seasons, implying that they had a significant shift when soil transformed from frozen to unfrozen. The season-driven shift in fungal communities may be partly explained by season-related changes in environmental factors (e.g., temperature and moisture) and biochemical variables (e.g., soil N availability and enzyme activity). Based on our findings, the intensified and extended soil frost associated with winter climate change might profoundly alter the phenology of soil fungal community in subalpine forests experiencing significant snowfall decrease.
Additional file 1: Figure S1. Seasonal dynamics of air temperature (a) (2 m above the ground surface) and soil temperature (5 cm depth) from November 2015 to November 2016 and soil moisture (b) in control and snow exclusion plots. The dot indicates snow depth during the winter. The asterisk indicates the sampling period. The small chart shows air and soil temperatures in winter. Figure S2. One-way ANOVA test bar plot for fungal phyla in (a) control (C) plots and (b) snow exclusion (SE) plots. DSP: deep snow period; ETP: early thawing period; MGS: middle in the growing season. * indicated a significant difference among sampling dates for a same snow treatment. Figure S3. Wilcoxon rank-sum test bar plot for fungal phyla in the deep snow period (DSP), the early thawing period (ETP), and the middle in the growing season (MGS). C: Control; SE: snow exclusion. * indicated a significant difference between snow regimes on a specific sampling date. Figure S4. One-way ANOVA test bar plot for fungal genus in (a) control (C) plots and (b) snow exclusion (SE) plots. DSP: deep snow period; ETP: early thawing period; MGS: middle in the growing season. * indicated a significant difference among sampling dates for a same snow treatment. Figure S5. Wilcoxon rank-sum test bar plot for fungal genus in the deep snow period (DSP), the early thawing period (ETP), and the middle in the growing season (MGS). C: Control; SE: snow exclusion. * indicated a significant difference between snow regimes on a specific sampling date.
Additional file 2: Table S1. Sequencing results and number of observed and estimated OTUs at the species level (70% 18S identity). Table S2. Total and relative number of sequences and operational taxonomic units (OTUs) distributed to different fungal orders.