Warmer winters result in reshaping of the European beech forest soil microbiome (bacteria, archaea and fungi) - With potential implications for ecosystem functioning

In temperate regions, climate warming alters temperature and precipitation regimes. During winter, a decline in insulating snow cover changes the soil environment, where especially frost exposure can have severe implications for soil microorganisms and subsequently for soil nutrient dynamics. Here, we investigated winter climate change responses in European beech forests soil microbiome. Nine study sites with each three treatments (snow exclusion, insolation, and ambient) were investigated. Long-term adaptation to average climate was explored by comparing across sites. Triplicated treatment plots were used to evaluate short-term (one single winter) responses. Community profiles of bacteria, archaea and fungi were created using amplicon sequencing. Correlations between the microbiome, vegetation and soil physicochemical properties were found. We identify core members of the forest-microbiome and link them to key processes, for example, mycorrhizal symbiont and specialized beech wood degraders (fungi) and nitrogen cycling (bacteria, archaea). For bacteria, the shift of the micro-biome composition due to short-term soil temperature manipulations in winter was similar to the community differences observed between long-term relatively cold to warm conditions. The results suggest a strong link between the changes in the microbiomes and changes in environmental processes, for example, nitrogen dynamics, driven by variations in winter climate.


INTRODUCTION
The global climate is changing rapidly and the global warming of 1.2-4.8C as well as a shift in precipitation regimes projected for this century will likely expose ecosystems to increased climatic stress (IPCC, 2013).Notably, climatic warming is projected to be stronger over land masses than over the sea and to increase towards higher latitudes (IPCC, 2013).In this context, most climate-change studies in northern temperate forest ecosystems have focused on the impact of increased drought and heat exposure on forest health (Allen et al., 2010;Lindner et al., 2010;Schuldt et al., 2020).However, the soil climate and the microclimate in the forest interior are generally decoupled from the macroclimate to a varying degree across seasons and forest types; therefore, it is not straightforward to make predictions for the fate of forest ecosystems purely based on average climate warming rates (De Frenne et al., 2021).In northern ecosystems, this decoupling can be very pronounced during winter in situations where a sufficient snow layer covers the earth, thereby insulating and moderating the soil temperature, that is, protecting it from frost (Kreyling, 2020;Studd et al., 2021;Sturm et al., 1997).However, climate warming is predicted to decrease snowfall in temperate ecosystems, which could paradoxically expose soil to increased frost frequency and severity (Groffman et al., 2001;Kreyling, 2010;Studd et al., 2021).Until recently, the winter season was considered a dormant period with little significance for ecosystems, and although still under-represented in climate change research, the winter climate is now receiving increased attention in temperate forest research (Kreyling, 2010(Kreyling, , 2020;;Studd et al., 2021).A number of dendroecological studies across northern temperate forest ecosystems have demonstrated a sensitive response in tree growth to mid-winter temperature variability, likewise large-scale forest experiments have indicated that coldinduced growth declines are related to soil cooling (Harvey et al., 2020;Pederson et al., 2004;Reinmann et al., 2019;Weigel et al., 2018Weigel et al., , 2021)).Finally, soil cooling has also been found to induce shifts in the composition of the forest understory vegetation (Kreyling et al., 2012).
Plant responses to colder soils have previously been explained by the direct effects of frost damage to roots resulting in a reduced nutrient uptake capacity (Comerford et al., 2013;Reinmann & Templer, 2016;Repo et al., 2021;Sanders-DeMott et al., 2018;Tierney et al., 2003) and subsequently decreased rates of decomposition of organic matter and increased nutrient leaching (Campbell et al., 2014;Joseph & Henry, 2008;Sanders-DeMott et al., 2018;Weigel et al., 2021).However, soil microorganisms are directly responsible for nearly all bio-geochemical transformations and thus play a key role in nutrient dynamics (Kallenbach et al., 2016;Nielsen et al., 2011).Thus, in the search for underlying causes for these reported forest ecosystem responses, one needs to consider the responses of the soil microbial community to alterations in temperature regime, in particular, increased temperature fluctuation, drought and frost exposure, which are considered severe environmental stresses for microorganisms (Schimel et al., 2004;Sorensen et al., 2016Sorensen et al., , 2018)).
The present study is conducted within a large study framework focusing on the distribution of European beech from central to cold-marginal populations (see section Experimental Procedures and Henry et al., 2021;Weigel et al., 2021).In brief, previous studies have documented that a core group of plants characterizes the forest understory at all investigated sites and only a small species turnover was seen across the distribution gradient (Weigel et al., 2019).Furthermore, differences in the abundance of these understory plant species were primarily explained by soil properties (organic matter content, pH and C/N ratio) and second by the average length of the growing season (1370-2140 C*h growing degree days above the threshold of +5 C) as well as moisture availability (mean annual precipitation 540-620 mm; Weigel et al., 2019).No effect was seen in the forest understory as a result of an experimental soil temperature manipulation in winter (See section Experimental Procedures and Weigel et al., 2021).In contrast, beech tree growth was shown to be sensitive to the experimental soil temperature treatments (where short-term soil cooling in a single winter was found to correlate with low tree stem growth in the following growing season; Weigel et al., 2018Weigel et al., , 2019)).Likewise, the short-term soil temperature manipulations changed soil processes at the study sites in the form of decreased litter decomposition rates, high over-winter nitrogen retention, and higher soil ammonium availability in the plots which had the lowest soil temperatures (Weigel et al., 2021).
Cold and sub-zero temperatures influence microbial activity by on one hand, simply slowing down metabolic rates and decreasing the rate of substrate processing in the soil (Schimel et al., 2004).On the other hand, intense freeze-thawing is a lethal threat to microorganisms due to frost dehydration and cell rupture (Jefferies et al., 2010;Schimel et al., 2007;Shchepin et al., 2014).A recent mesocosm experiment found the strong response to altered winter temperature conditions in the soil microbial community composition, which was interlinked to changes in the nutrient dynamic and plant growth in the following season of two common temperate plant species (Calluna vulgaris and Holcus lanatus; Borg Dahl et al., 2021).Additionally, plant-soil interaction models predict a faster response to environmental changes among soil microbes than among primary producers, for example, forest trees, with the former likely undergoing a community reordering prior to the latter (Bardgett et al., 2013).Documenting microbiome community changes has proven difficult due to the high spatial and temporal complexity of soil microbial communities and the physicochemical heterogeneity of soil.With the additional difficulties related to winter sampling, many previous studies have been constrained to low resolution, for example, monitoring gross soil microbial activity or total biomass (Dur an et al., 2014;Sorensen et al., 2016Sorensen et al., , 2018)), or such studies only compare between singlepoint microbial community compositions in summer versus winter (Monson et al., 2006).With modern molecular technologies, high-resolution microbial studies are possible, but there has hitherto been a strong focus on high-organic boreal and arctic soils (including permafrost conditions) due to the carbon storage found in these soils (Crowther et al., 2016).In addition, many studies have investigated the short-term response to elevated temperatures using laboratory incubations ( Čapek et al., 2015;Elberling et al., 2013;Schostag et al., 2019;Tveit et al., 2015).In contrast, limited knowledge exists on how winter climate change scenarios in temperate forest ecosystems, that is, under rather moderate temperature changes, might affect the soil microbial community in terms of species composition and abundance.Urakawa et al. (2014) transplanted soils among temperate forest sites, which either naturally experienced non-frozen or frost-effected winters in Japan.The authors found the nitrogen dynamics (nitrogen mineralization, nitrification and denitrification) highly affected by winter climate (with legacy effects into the following growing season) and that both the magnitude, that is, the absolute sub-zero temperature, and frequency of freeze-thaw cycles were important explanatory parameters of their results.These findings indicate effects on the functional microbial groups, such as ammonia-oxidizing archaea and bacteria (Leininger et al., 2006), as well as denitrifiers (Philippot et al., 2009).In line with the findings of Urakawa et al. ( 2014), Kreyling et al. (2020) found in a laboratory experiment that both increased magnitude and frequency of freeze-thaw cycles led to increased nitrate, ammonium, and phosphate release almost exclusively in soils originating from snow rich (i.e. less soil frost affected) sites within the present study frame.Thus, here we investigated the bacterial, archaea and fungal soil community compositions in all experimental plots within the study framework of the European beech distribution.Apart from the wider aim to give a first description of the beech forest microbiome and the expectation to identify a so-called 'core microbiome', that is, microbial taxa common in microbial assemblages associated with a certain habitat, here the northern European beech forest (Neu et al., 2021;Shade & Handelsman, 2012), we specifically hypothesize (1) that the shift of the soil microbial community from central to cold-marginal beech populations resembles the floristic gradient observed for the forest understory plant community, with a gradual change in microbiome composition instead of a complete community turnover across the study area.Further, we hypothesize (2) that winter soil temperatures have a role in shaping the microbial community composition and (3) that especially shortterm soil temperature variation in winter influences the microbial community composition more than previously observed for the forest understory vegetation.Finally, we investigate whether such short-term changes can be linked to previously observed ecosystem responses, for example, nitrogen availability.

Site characteristics
The study sites were established on an east-west winter temperature gradient across 500 km (ΔT = 4 C decrease towards the east) in nine temperate deciduous forests in North Germany and North Poland.The managed (selective thinning every 5-10 years), semi-natural to natural forests are mono-dominated by mature European beech trees (Fagus sylvatica L.) and range from the central to the cold-marginal populations within the beech distribution (Figure 1A,B).We selected forests with a similar hall-like forest structure and similar tree stature with mature trees ranging 37-56 cm in diameter, 27-40 m in height, and 76-167 years of age on average across sites (Weigel et al., 2018), which represents the intermediate state of typical beech forest succession in the study area (Leuschner & Ellenberg, 2017).The last thinning at the sites was ≥5 years prior to the experiment and no forest operation took place during the experiment and the experimental plots were spaced one tree length apart from any canopy gaps.The climate in the study area is cool to cold temperate with 530-630 mm of mean annual precipitation and 5.6-8.8C mean annual temperature ('Clima-teEU' v4.63 data;Hamann et al., 2013;Wang et al., 2021).Snow can occur from Novemberto March (April) at all sites and the winters become markedly snowier towards the east (Weigel et al., 2021).The soil at all sites is of the type sandy Cambisol, which formed on glacial deposits in the Pleistocene lowlands.The sites differed slightly in the clay, silt and sand content (measured in laser particle sizer, Analysette 22 microtec, Fritsch, Idar-Oberstein, Germany), with DE, NZ and BH (see Table 1 for an explanation on the site IDs) being the sandiest soils (and with the lowest amount of clay).For the latter two sites (NZ and BH), the organic matter (OM; measured as the weight fraction loss during combustion at 550 C for 12 h) and the soils C:N ratio (measured with vario EL III Element Analyzer, elementar Analysensysteme GmbH, Hanau, Germany) were among the lowest, while the DE site had the highest C:N ratio and high OM content, making DE the most unique site (Table 1).A weak pH (measured in a 0.01 M CaCl 2 solution in a 1:5 soil: solution ratio 12 h after shaking) gradient existed across the sites, with DE having the lowest pH (3.0) and BB the highest (pH 3.5).The soil temperature was monitored at 5 cm depth (between the O [litter layer] and A [topsoil] horizon) every 30 min ('TMC20-HD" temperature sensor connected to a 'HOBO UX120-006M Analog Data Logger'; Onset Computer Corporation, Bourne, USA).Based on the average and minimum soil temperature, the sites rank as following (coldest to warmest); KO, GR, DE, WE, KA, BB, NZ, HH and BH with slightly ambiguous positions in the middle-ranked sites but significant differences between BH and KO, this order will be kept throughout the article as well as the corresponding colour code as shown in Figure 1C-F.

Experimental set-up and soil sampling
At each site of this gradient-design study (replication across many sites, instead of many replications at few sites, Kreyling et al., 2018) three fenced experimental plots were established in autumn 2016.Each plot was a circle with a 10 m diameter, corresponding to a plot size of approximately 80 m 2 , and with a representative mature tree in the center.During the treatment period (November 2016-March 2017), one of the three plots at each site was assigned for the snow exclusion treatment (SÀ), where snow was blocked from falling to the ground by a transparent plastic cover (0.2 mmthick, UV-transmissive 'Lumisol clear'; Hermann Meyer KG, Rellingen, Germany) that was carried in 1.8 m height by a wooden roof construction.The plastic cover was perforated and thus permeable for snowmelt water F I G U R E 1 Treatment effects on soil temperature characteristics.Location of the study area in Europe (A) with the study sites situated along a strong west-east winter temperature gradient (B).Measured variables (C-F) from the experimental period 2016/2017: 1 December to 12 March for each site (see colour legend and see Table 1 for an explanation on the site IDs) and each treatment in the order SÀ, C, S+.Significant differences between treatments and control plots (SÀ vs. C, S+ vs. C for each variable) across all sites are indicated with ***p < 0.001 (paired t-test), ns, non-significant.Ave.minimum soil temperature refers to the average of all daily minimum temperatures during the study period.and winter rainfall.One plot was assigned as an insolation treatment (S+), where additional insulating snow cover was simulated by a gardening fleece cover (70 g/ m; Hermann Meyer KG, Rellingen, Germany).One plot was assigned as a reference for ambient conditions (control 'C').At the site with on average coldest winter air temperatures (Kartuzy, KA), the treatments were triplicated (n = 9) and these additional plots were used for the analysis of short-term treatment responses.Thus, in total n = 33 treatment plots were established along the gradient.Directly after the end of the treatment period, soil samples (ca.20 g fresh soil) were collected in a time window of 10 days in March 2017 from the topsoil in three evenly dispersed (northeast, southwest-south, and northwest in 3 m distance from the plot center) replicates per plot, corresponding to n = 99 individual samples in total.The samples were immediately placed on ice and stored at À20 C until further processing.
Due to the SÀ treatment only resulting in marked soil cooling at the KA site (see Results section 'Soil temperature responses to the snow manipulation experiment'), the following analysis will distinguish between, first, the short-term effects of winter soil temperature variability in a subset of our data.We did so by comparing only the single-winter SÀ and the S+ treatments at KA (n = 6 plots).Second, in a different subset, we compare the ambient condition C plots between the site with the warmest (site BH) and coldest (site GR) soils on the gradient, which should represent the long-term effects of soil microbial acclimatization to natural average soil temperature conditions.GR likely experiences a more sporadic snow cover and thus a colder and more fluctuating-and frost-dominated soil temperature (as was the case in the experimental year, Figure 1C-F), although the sites KA and KO rank as colder than GR when considering air temperature.In addition, a high similarity in soil physiochemical parameters exists between BH and GR (Table 1), which is crucial for an evaluation restricted to single-site comparisons.

DNA extraction and library preparation
Using a micro scale, 0.25 g of each soil sample was added to a PowerBead tube (QIAGEN, Valencia, CA, USA) and placed into a FastPrep homogenizer (MP Biomedicals, Irvine, CA, USA) at 5 m s À1 for 45 s.The DNA extraction was done using QIAGEN DNeasy PowerSoil Kit following the manufacturer's protocol.DNA content was measured using DeNovix DS-11 FX Spectrophotometer and ranged from 50 to 100 ng DNA/μL.For PCR reactions, DNA extracts were diluted to a concentration of ca. 10 ng μL À1 and 1 μL was added to the PCR reaction using GoTaq (Promega, Madison, WI, USA) buffer and polymerase (final volume of 25 μL).

Bioinformatic processing: Taxonomic and functional classification
Using R v.4.1.2sequence reads were quality filtered and amplicon sequence variant (ASV) tables were generated via the Dada2 pipeline (Callahan et al., 2016) and available form suppl. S1.For the fungal samples, the DNA read quality of the forward (R1) reads was poor, and overlapping of the forward and reverse reads was not possible.Therefore, an ASV table was generated based only on the reverse (R2) reads, which were found to be of good quality up to 220 bp (suppl.S2).For bacteria/archaea, the forward and reverse reads were combined using FLASH (Magoč & Salzberg, 2011) prior to quality filtering.In brief, the following quality filters were applied for 16S and ITS: 16S (overlapped reads): maxEE = 2, minLength = 50, trim-Right = 20, trimLeft = 19 (additional parameters were run as default), ITS (R2 reads): truncLen = 220, maxEE = 2, trimRight = 20 (additional parameters were run as default); 85 and 90 samples successfully passed the quality control for fungi and bacteria/ archaea respectively (for detailed documentation see suppl.S3 and S5).The three failed samples for 16S were; one failed to amplify (not processed), one had a very low richness (241 ASVs vs. 1014 ± 283 for the remaining samples) and one was strongly dominated by two ASVs with very high abundance (26%) obscuring the community structure in this sample.For fungi, seven samples had <1000 reads and were removed (on average 16,858 ± 5576 reads were found for the remaining samples) and one had low ASV richness (42 ASVs vs. 123 ± 34 for the remaining samples).Negative controls were found satisfying with <500 reads for 16S (compared to 67,880 ± 15,662 reads per sample) and <20 reads for ITS.ASVs were assigned taxa using Usearch v. 11.0.667(Edgar, 2016) with the RDP 16 S v16 database (Cole et al., 2014) and the UNITE database v. 7.2 (Abarenkov et al., 2010) as reference for bacterial/archaea (16S) and fungal (ITS) reads, respectively.For 16S only ASVs assigned to the domains 'Bacteria' or 'Archaea' were kept for the downstream analysis and within the Bacteria ASVs assigned to 'Chloroplast' were removed.For both datasets, ASV counts were normalized using CSS (Cumulative Sum Scaling; Paulson et al., 2013).The ASV tables and taxonomic assignments are available in suppl.S1.The FUNGuild database (v.2) was used for functional classification of the fungal community (Nguyen et al., 2016).Only ASVs assigned to family level or deeper could be classified.Only classifications with a high confidence score ('highly probable' or 'probable', i.e. omitting 'possible') were used.The following ecological groups were found: mycorrhizal fungi (arbuscular, ecto-, ericoid and orchid mycorrhizal), endophyte, lichenized, plant pathogen, saprotroph, fungal parasite.For ambiguous classifications the following ranking was applied: potential ectomycorrhizal (EcM), potential wood saprotroph, potential fungi/lichen parasites, potential plant pathogen and finally 'animal-associated fungi', in most cases, the alternative function was 'saprotrophic' or for the pathogens 'endophytic'.
Likewise, FAPROTRAX (1.2.4 with python 3.7) together with provided script and database (Louca lab., Department of Biology, University of Oregon, Eugene, USA) was used for functional profiling of the bacterial/ archaea community.As with the fungi only ASVs assigned to family level or deeper could be classified.In addition, ASVs belonging to the archaeal phyla thaumarchaeota were further assigned to (AOA) amoA clades against a database containing 16S rRNA gene sequences associated with their amoA gene counterparts, according to Wang et al. (2021).This enables access to a novel and highly resolved amoA taxonomy (Alves et al., 2018).The assignment was done with UCLUST based on an 85% similarity.

Data analysis and statistics
Basic statistical analyses were performed in R v. 4.0.3(R Core Team, 2020) using mainly the packages 'vegan' (Oksanen et al., 2020), and 'phyloseq' (McMurdie & Holmes, 2013).'ggplot2' was used for visualization (Wickham, 2016).Environmental treatment effects are all reported as averages with ±standard deviation, and statistical significance differences (p < 0.05) between means were tested using paired t-test.Distance-based redundancy analysis (dbRDA) was carried out using Bray-Curtis dissimilarity matrix.Investigated variables were 'Site (KO, GR, DE, WE, KA, BB, NZ, HH and BH)', 'treatment (C, SÀ and S+)', 'Shannon Diversity index for the vegetation', 'sample score for 1 st NMDS axis for the vegetation composition' and 'sample score for second NMDS axis for the vegetation composition' (obtained from Weigel et al., 2021), 'soil frost-degree-hours', 'average soil temperature', 'minimum soil temperature', 'freezethaw-cycles', 'NO3-N availability', 'NH4-N availability' (measured as absorption to resin plates throughout the experimental winter, for details see Weigel et al., 2021), 'Green tea decomposition', 'Tree stem growth'.The partitioning of variance on the significant axis was reported (p < 0.05, PERMANOVA with 999 permutations).For the dbRDA biplot ordination, the sample scores (weighted sums of species/ASVs) and the correlation biplot scores of the environmental parameters were used for plotting the ordination (in which angles between all vectors reflect linear correlation).The contribution of individual variables was tested with PER-MANOVA adding variables sequentially in the order listed above.Additional calculations were made for the site alone and site-treatment interactions as explanatory variable.The similarities between the fungal and bacterial communities between samples were tested with Mantel tests of a Pearson's correlation of Bray-Curtis dissimilarity for bacteria/archaea and fungi as well as with a Procrustes test on two-dimensional NMDS ordinations of each dataset with the function protest() with 999 permutations in the R package 'vegan' (only shared samples included, n = 83).Heat trees were made using the R package 'metacoder' (Foster et al., 2017) showing only taxa with a summed abundance >0.005 (bacteria/archaea) and >0.05 (fungi), colour gradient, node size and -label corresponding to relative abundance.Core members of the microbial community are defined as ASVs occurring in a minimum of one sample at all nine sites.
To compare the community response to long-term temperature conditions, the control (C) plots for warm (BH) and cold (GR, in terms of soil temperature) were compared (n = 6 samples) by calculating the log2-fold change for ASVs, summarizing the abundance from genus to phyla level.Only taxa with a presence in >50% of the samples (i.e. two of three or four of seven samples for long term or short term, respectively) and an average abundance >0.5% per treatment were reported.Differential abundance between warm-and cold conditions was tested with the Student's t-test.The analysis script is available in Suppl.4.

Soil temperature responses to the snow manipulation experiment
While the insulation in the S+ treatment significantly increased average-and minimum soil temperature and decreased the frost degree hours across the gradient (Figure 1C,D), the SÀ treatment had a little-to-no effect (and in some cases, even a reverse mini-greenhouse effect was observed) at most sites.The winter was unusually warm and snowfall was lacking across the gradient in 2016/2017 (Figure 1 and Weigel et al., 2021).Thus, the SÀ treatment did not result in a uniform temperature treatment across sites.Only at the KA site was an insulating snowpack present for periods of time across the winter, and here significant (p < 0.05) differences between treatments were seen for all variables (Figure 1C-F, dark grey).

Fungal and prokaryotic core microbiomes of beech forest soils
In total, 27,392 and 2990 ASVs were obtained for bacteria/archaea and fungi, respectively (suppl.S6).Although individual ASV count largely differed, the number of unique taxa was in a similar range between bacteria/archaea and fungi, for example, 243 and 200 genera were represented in each community.For fungi 12 of the existing 16 phyla (Wijayawardene et al., 2018) were identified, for bacteria/archaea 26 phyla were represented of the currently 41 phyla recognized word wide (https://lpsn.dsmz.de/).
For both fungi and bacteria/archaea, some ASVs were found in samples at all sites (on average occurring in 50 ± 20/54 ± 18 samples for 16S/ITS, respectively) corresponding to ca. 1.5% of the ASVs (434/27,392 ASVs and 44/2990 ASVs for bacteria/ archaea and fungi, respectively; Figure 2A,B), however accounting for 59.6% ± 0.1% and 30.7% ± 0.8% of the total bacterial/archaeal and fungal abundance, respectively.These ASVs were considered core members of the beech forest microbiome and their taxonomic composition along with the overall composition is shown in Figure 3A-D.

Impact of site and treatments on fungal and prokaryotic microbiomes
The bacterial/archaeal and fungal community composition showed a similar assembly between samples (Pearson's product correlation of 0.43, Mantel statistic p < 0.001 and Procrustes m12 squared = 0.5, R = 0.7, p = 0.001, suppl.S5) and the variance explained by the investigated variables was similar for the two data sets (total constrained variance through explanatory variables was 39.1% for bacteria/archaea and 37.2% for fungi), with 'site' explaining most of the variance in the communities' compositions for both (27.1% for bacteria/archaea and 23.3% for fungi PERMANOVA, p < 0.001; suppl.S5).The experimental treatments alone could only explain a minor part of the community variance and were only significant for fungi (2.6%, p < 0.05).However, the factor 'treatment' appeared to have different effects across the sites (seen as a strong interaction effect between treatment and site, 18.8% for bacteria/archaea and 19.3% explained variance for fungi p < 0.001; suppl.S5).This corresponds to the little-to-no effect of the cold (SÀ) manipulation on soil temperature at most sites due to unusually warm winter climate in the experimental year (see above).
The first five RDA ordination axes could significantly constrain variance in both communities, with the first two axes explaining 14.9% and 12.3% of the variance for bacteria/archaea and fungi, respectively (Figure 4A,B and suppl.S5).For both bacteria/archaea and fungi, vegetation properties of the forest understory (Shannon Diversity index and NMDS axis 1 and 2 of the plant community composition) were strongly related to the community profiles, that is, strongly correlated with the first two dbRDA axes (see biplot scores for constraining variables, suppl.S5).The warmer sites (BB, BH and NZ with higher average and minimum soil temperature) were associated with increased availability of NO 3 -N (which for fungi was strongly associated with the first dbRDA axis, suppl.S5), whereas the relatively colder sites (KO, KA and WE) were associated with increased availability of NH 4 -N and the vegetation NMDS-axis 2. The frostdegree-hours and the number of freeze-thaw cycles showed only weak associations to both community compositions.
The site DE showed especially strong separation from the remaining sites, indicating a unique composition of both fungi and bacteria at this site, corresponding to the most unique soil properties among the study sites (highest CN ratio and sand content, Table 1).For both communities (bacterial/archaeal and fungal) sam-F I G U R E 2 Overview of amplicon sequence variants (ASVs) distribution.First (inner) ellipse represents the average number of ASVs (± S.D.) found in replicates within a treatment plot, second ellipse: ASVs found in replicated plots within the KA site, third ellipse: ASVs found in treatment plots (C, SÀ and S+) per site.S, shared numbers (found in minimum two out of three) replicates/plots/treatments; T, total numbers (found in any one of three replicates/plots/treatments).Fourth ellipse: ASVs found in all and five out of nine sites.(A) Bacteria/archaea.(B) Fungi.
ples from the site KO varied a lot among each other.Samples from the replicated short-term treatment design at site KA tended to form clusters at positions corresponding to the overall cold-warm gradient orientation, that is, relatively warm treatment samples clustered at dbRDA1 > 1 (like the warmer sites), while cold and ambient treatment clustered towards dbRDA1 < À1 or dbRDA2 > 0 (like the colder sites, Figure 4), this was only seen for the bacterial community (i.e.no treatment specific clusters were seen for fungi).

Short-term versus long-term temperature response of microbiomes
Several bacterial taxa had a significantly different abundance between warm and cold long-term conditions (Figure 5A), while many taxa had a similar distribution under the short-term conditions (log2 values at zeroline, p > 0.1).The bacterial taxa which did show a differentiated abundance under short-term conditions (although the test showed only a trend for most of them; p < 0.1; Figure 5A) often showed similar (significant) patterns under long-term conditions, too.These were members of the Bradyrhizobiaceae (Alphaproteob acteria) and Deltaproteobacteria which both were found in higher relative abundance under warmer conditions, the pattern persisted at lower taxonomic levels where the order Myxococcales (and family Polyangiaceae), as well as the genus Bradyrhizobium, had a similar response to soil temperature.Thus, in accordance with the prior ordination (dbRDA) analysis, the abundance pattern between long-and short-term conditions showed some level of congruency.
For fungi, few significant responses were found and more contrasting abundance patterns were seen between long-and short-term conditions (Figure 5B).In the short term, a significant response of the dominant Basidiomycota order Agaricomycetes was seen (higher relative abundance under warmer conditions).Among observed trends (p < 0.1) were a lower relative abundance of the genus Xerocomellus (Basidiomycota, Boletales) under cold conditions, but a higher relative abundance of the genera Umbelopsis (Mucoromycota) and Mortierella (Mucoromycota).

Bacteria/archaea
Of the 2584 bacterial/archaea ASVs assigned to family level or deeper 1635 were functionally classified via FAPROTAX (an additional 373 ASVs of higher taxonomic levels were also assigned to a function, this was mainly members of the phyla Chlamydiae, which were assigned to the functional group 'intracellular parasites').In total, 45 unique functions were found (suppl.S1).Among the classified ASVs, each ASV was assigned to a minimum of one function and a maximum of 13 functions (two ASVs) and the median was two functions per ASV.Assigned ASVs accounted for 20.3% ± 3.6% of the relative abundance per sample.By far the most dominant classification was 'Chemoheterotrophy' (an-and aerobic) with 1346 ASVs accounting for 7.2 ± 1.5% relative abundance per sample, followed by 'intracellular parasites' (635 ASVs; 1.1 ± 0.4%), 'fermentation' (177 ASVs; 1.5 ± 0.4%) and 'nitrate reduction' (167 ASVs; 1.4 ± 0.4%).Of the core bacterial community (434 ASVs), only 47 ASVs could be assigned to a function.The core functional profile resembled the overall functional profile with mainly 'Chemoheterotrophy' (43 ASVs), 'fermentation' (13 ASVs) and 'nitrate reduction' (13 ASVs).Four additional functions were identified, each represented by a F I G U R E 4 Sample similarity in community composition and main environmental parameters.Shown are the first and second axis of the distance-based RDA ordinations for the bacteria/archaea (A) and fungi (B).Dots are coloured by site and shaped by condition (see legends) and scaled to the number of amplicon sequence variants (ASVs) per sample.Environmental variables are plotted as overlaying vectors (angles between vectors reflect linear correlation); FTC, freeze-thaw-cycles; FDH, frost degree hours; NH 4 , bulk soil ammonium availability; NO 3 , bulk soil nitrate availability; Shannon Div., diversity index for the plant community.Tave, average soil temperature, Tmin, minimum soil temperature.Vegetation NMDS1 and 2 are axis scores for the plant community (see Experimental Procedures for details).single core community ASV, which were: 'aerobic ammonia oxidation' (taxonomic identity of ASV: Archaea; Nitrososphaerales), 'nitrification' (Archaea; Nitrososphaerales), 'nitrogen fixation' (Alphaproteo bacteria; Bradyrhizobium) and 'dark hydrogen oxidation' (Alphaproteobacteria; Labrys).The taxonomic F I G U R E 5 Response to relative warm and cold soil temperatures under long-term (cold GR, warm BH) and short-term (KA site; S+ and SÀ) conditions of bacteria/archaea (A) and fungi (B).Taxa from each taxonomic level (sorted from phyla to genera) with a presence in >50% of samples (two of three or four of seven samples for long-and short-term, respectively) and an average relative abundance >1.0%.Dots are scaled to average relative abundance (see scale bar) and coloured by significance (see legend).
clade classification of the Thaumarchaeota ASVs (AOA) showed that the NS-gamma clade was by far the most abundant, followed by NT-alpha, NS-alpha and NS-Zeta (following Wang et al., 2021, suppl. S1).

Under long-term conditions
'Intracellular parasites' were less abundant under warmer conditions (p = 0.03; Suppl.S5).In addition, the functional classifications 'photoheterotrophy', 'phototrophy' and 'predatory or exoparasitic' were only found under warmer conditions in the long-term comparison (all with <0.1% average abundance).For the latter ('predatory or exoparasitic') the functional classification exclusively consisted of ASVs assigned to families (or genera) within the orders Myxococcales and Bdellovibrionales (Deltaproteobacteria).ASVs assigned to one of these two taxonomic orders were also found under long-term cold conditions (see Figure 5A), however not assigned to family (or deeper) level and thus not functionally classified.

Fungi
Of the 1169 fungal ASVs assigned to family level or deeper, 918 were functionally classified via FUNGuild.These ASVs accounted for 59% ± 14% of the abundance per sample.The most dominant groups were EcM fungi (ca.40% of the total community), free-living saprotrophs (ca.10%), ericoid mycorrhizal (ca.5%) and wood degraders (ca.3%, suppl.S5).None of the functional groups (with an average abundance >1%) showed a difference in abundance between relative cold and warm conditions.

DISCUSSION
The core microbiome of beech forest soil We investigated the bacterial/archaeal and fungal soil community composition within the north-eastern range of European beech forest ecosystems across a winter climate gradient.For both communities, we identified a number of taxa present at all sites and considered them to represent core members of the beech forest soil microbiome.The taxonomic identity of these taxa was in accordance with findings from other studies of the beech forest soil microbiome (Choma et al., 2020;Jeanbille et al., 2016;Tedersoo et al., 2014).From the functional classification of the prokaryota, well-known nitrite-oxidizing bacteria (NOB) and ammonia-oxidizing archaea (AOA) were identified, namely the archaeal class Nitrososphaeria (Thaumarchaeota) and the bacterial genus Nitrospira.The AOA was part of the core microbiome and, accordingly, studies have found AOA to be the predominant ammonia-oxidizer in various soils (Leininger et al., 2006), especially under low ammonia availability (Aislabie et al., 2013) and low-pH conditions (Prosser & Nicol, 2012;Zhang et al., 2012).Likely, the strongly acidic soil conditions (pH 3.0-3.5) in this study were selected for AOA as the primary ammonia oxidizers.In accordance with this, the NS-gamma and NT-alpha clades predominated the AOA community, and these two AOA clades rather than other AOA clades are dominant in acidic soils and known to persist under low ammonium concentrations (Alves et al., 2018;Wang et al., 2021).These findings support the identification of AOA as key members of the forest soil microbiome.However, much more dominant (and also a core member of the community) were members of the genus Opitutus (Verrucomicrobia), from the functional classification these were also assigned to taxa carrying out nitrification although the specific pathway for these taxa remains unknown (Baskaran et al., 2020).Also, taxa carrying out denitrification were identified among the core microbiome, namely the genus Schlesneria (Planctomycetes), for which only one species is currently described (S.paludicola isolated from acidic peat bogs and capable of utilizing nitrate; Kulichevskaya et al., 2007), and Paracoccus (Alphaproteobacteria), for which many described species exist, for example, the highly studied denitrifier P. denitrificans (the amplicon length of the ASV does not allow for a species identification in this current study).
Among the fungi, the identified core members were mainly EcM fungal (e.g.Genea hispidula), saphrotrophs (Mortierella angusta and Umbelopsis dimorpha) and specifically efficient beech wood degraders (e.g.Phallus impudicus, Heilmann-Clausen, 2001).In addition, many of the identified EcM fungi have been reported as specific symbionts with European beech (e.g., Genea hispidula, Lactarius subdulcis, Byssocorticium atrovirens and Clavulina sp., Russula sp., Tomentella sp. and Cortinarius sp.; Buée et al., 2005;Köhler et al., 2018).This is not surprising as EcM fungi play a key role in the forest ecosystem, where they account for up to 50% of the microbial biomass, occurring in the form of dense hyphal mats often with a low pH (Baldrian, 2016;Kluber et al., 2010).Bertaux et al. (2005) found EcM fungi to host endosymbiont bacteria of mainly alphaproteobacteria, a large and highly diverse class of Bacteria, but the authors were not able to identify the endosymbionts to a deeper taxonomic resolution.In the present study, we find alphaproteobacteria among the most dominant bacterial class in the community (and identified as a core member of the forest microbiome).Among these, we found the order Sphingomonadales (mainly the genus Sphingomonas) most abundant, which has been associated with EcM fungi in other studies (Zimmer et al., 2009) and thus could be a likely candidate as the EcM endophyte described above.Among other identified, alphaproteobacterial were the well-known N-fixing plant-symbionts Rhizobiales (incl.the genera Bradyrhizobium, Rhizobium and Methylocella, the latter, a facultative methanotroph found in acidic forest ecosystems, Dunfield et al., 2003).

Strong interlinkage between microbiome, vegetation and soil physicochemistry
A major part of microbial community variation in our study was due to differences across sites.Despite the large range of the study gradient, the differences in microbial community composition were not due to the geographical distance between sites per see, but rather due to the edaphic and floristic (previously reported, Weigel et al., 2019 and Table 1) differences between them.Furthermore, a similar environmental response between bacteria/archaea and fungi was observed (correlation in community similarities), which is in line with the well-known strong link between plants and microbes, especially in systems dominated by symbiotic relationships as ours (Bahram et al., 2015;Berg et al., 2016) and their interlinked dependence on nutrient availability ( Čapek et al., 2018) as also indicated from the clear association between vegetation variance and the availability of especially NH 4 -N.Thus, in line with our first hypothesis, our results support prior observations, that the plant community composition is among the main drivers for bacterial/archaeal and fungal soil community assembly in a temperate (alpine) forest gradient (Borg Dahl et al., 2019).Likewise, Choma et al. (2020) found a strong and rapid change in the bacterial community composition in a beech forest as a result of an experimental change in pH of ca.0.5 unites (pH 4.3-3.7).The same authors found only minor changes in the fungal community in their experiment.In accordance, we observed the communities of both bacteria/archaea and fungi to be most distinct from each other between the sites with the lowest (DE, pH 3.0) and highest (BB, BH, pH 3.5) pH.
Microbiome response to warmer winter temperatures: Short-and long-term effects From the dbRDA ordination, a separation between relatively cold and warm sites was seen.Furthermore, the similarity between how the prokaryotic communities related to long-term temperature differences between sites and how communities responded to the short-term snow removal experiment suggests that the regulation of soil temperature by snow cover during the cold season truly has a part in shaping the soil microbial community in cool to cold temperate beech forests, as hypothesized (H2).Generally, we found more significant responses among the bacteria/archaea than fungi, which is in line with another finding of bacteria/archaea showing higher responsiveness than fungi to short-term environmental changes, as shown for beech forest soil microbiomes (Choma et al., 2020).Such short-term responsiveness is in line with our third hypothesis, namely that microorganisms might undergo changes at a different time scale than, for example, plants (no response was seen for the forest understory composition; Weigel et al., 2021).Such a decoupling in response might affect ecosystem functioning and the availability of plant nutrients (Bardgett et al., 2013;Yanai et al., 2004).In accordance, a reduced tree growth was observed from the short-term soil cooling experiment (SÀ affected plots) as well as from soil frost-affected sites along the gradient (Weigel et al., 2021).Among the bacteria, the core-member genus Bradyrhizobium had a consistently reduced relative abundance under short-and long-term cold conditions.Recent studies of this genus have shown that it is one of many free-living taxa that are dominant in forest soils (VanInsberghe et al., 2015) and that some are capable of N-fixation under non-symbiotic conditions (Tao et al., 2021), making them potentially important players in the nitrogen cycling of forests.
Another core community member of Bacteria which had a consistently significant lower relative abundance under short-and long-term cold conditions were members of the order Myxococcales (families Polyangiaceae and Myxococcaceae).These taxa (commonly referred to as myxobacteria) might be of particular relevance for the microbial community dynamics, due to their unique ecology and predatory behaviour.A recent study showed that myxobacteria were the predominant bacterial predator in mineral soils, with a higher abundance than both protozoa and nematodes (Petters et al., 2021), suggesting an important role in the microbial food web and thus in nutrient turnover.Our findings thus suggest alterations in the trophic interaction resulting from changes in the winter soil temperature regime.
We found no consistent responses to short-and long-term temperature variation in neither the taxonomic-and functional profile of the fungal community (only a slight decrease in the major phylum Basidiomycota as a result of the SÀ treatment), which partially contrasts our second hypothesis.Species assembly in ecosystems dominated by symbiont species such as the forest soil microbiome has been shown to not follow general distance-decay rules due to the selective impact of the plant-host community (here trees), acting as a buffer against changes in microclimate (Goldmann et al., 2016), which might be part of why we do not observe consistent patterns.Albeit statistically non-significant (or only near-significant, p < 0.1), some short-term responses were observed among the fungi, which might be explained by taxon-specific differences in lifestyle.Examples of these were the decrease in the relative abundance of the genera Xerocomellus, and the increase of Umbelopsis and Mortierella as a result of the SÀ treatment.Xerocomellus (exclusively X. pruinatus, one of the fungal core microbes) is an EcM fungi shown to grow primarily with vegetative (clonal) growth and form potentially large (>100 m 2 ) individuals connected to several host trees (Fiore-Donno & Martin, 2001;Šutara, 2008).The change in abundance likely is a result of a reduction in growth under the frost exposure, rather than increased growth in response to the slight increase in temperature under the S+ treatment, which is in accordance with the reduction of beech tree growth due to soil cooling in the experiment (Weigel et al., 2021).It is possible to imagine that frost has either (i) caused damage/death to the hyphae or (ii) resulted in a reallocation of mycelium growth, that is, in an attempt to avoid/escape the frost treatment.Similar patterns were observed for plant roots in the experiment, where the growth of fine roots shifted deeper into the soil following increased frost exposure (Weigel et al., 2021).Although fungi have been hypothesized to be more susceptible to frost damage than other microbes due to the predominant filamentous growth contra single-celled microbes (Schimel et al., 2007), EcM fungi of Betula and Alnus trees have been shown to withstand hard frost events (Kilpeläinen et al., 2016).However, Tibbett et al. (2002) found a poor recovery of EcM fungi after freezing events if the temperature prior to freezing was high, that is, if exposed to large temperature amplitude the frost-tolerance of EcM fungi decreased.Temperature fluctuations increasing in frequency and amplitude are expected for soils when the insulating snow cover declines, such as in our study where the lack of snow cover caused a more than 3-fold increase in freeze-thaw cycles and a more than 4-fold increase in frost-degree-hours at the coldest site (Figure 1).The documented association to multiple hosts in combination with the large individual sizes described above for X. pruinatus might be factors that help secure their survival (Lian et al., 2006;Šutara, 2008).Similarly, a high EcM biomass turnover as well as a small scale (0-20 cm) mobility in the course of a year have been documented for EcM fungi in an old-growth mixed-conifer forest (Izzo et al., 2005).Contrary, Umbelopsis and Mortierella (both previously classified under the former phylum Zygomycota), with the popular name 'sugar fungi' (utilizing mainly substrates rich in simple carbohydrates) are fast-growing 'moulds' that benefit from easily degradable substrates, such as recent microbial necromass (Richardson, 2009).Lindahl et al. (2010) found members of Mucoromycetes (former phylum Zygomycota) to increase in abundance upon the death of EcM fungi in an in situ root disruption experiment under conifers.A similar dynamic was observed between ericoid mycorrhizal fungi and members of Mucoromycota in a recent mesocosm experiment with the shrub Calluna vulgaris where the soil was exposed to winter pulse-freezing treatments, resulting in a strong decrease of the mycorrhizal symbiont and a high increase in Mucoromycota sp.(Borg Dahl et al., 2021).Taken together with our observations, this suggests that a similar dynamic occurs in our SÀ treatment; unfortunately, no quantitative measures of EcM biomass are available to confirm this hypothesis.
Both taxonomic and quantitative changes in the microbiome will affect the nutrient dynamics.In the functional profiles of bacteria/archaea, we found a slight increase in ASVs potentially capable of 'nitrification' along with increased NH 4 + availability after the short-term experimental snow removal and soil cooling of a single winter.Similarly, the group of potential denitrifying bacteria had a higher relative abundance under short-term warm conditions co-occurring with increased NO 3 À availability.It is reasonable to assume a connection between these observations, but a direct cause and effect cannot be concluded with the data at hand.In the future, a targeted approach could be applied to shed light on these soil processes, for example, by directly targeting process measures and gene expressions.
In summary, our results indicate a strong linkage between soil chemistry, plant and microbial community composition, with consequences for environmental processes and nutrient dynamics.In cold-marginal, previously snow-affected beech forests, warmer winters can be expected to cause a reorganization of the soil microbiome (especially bacteria) into a community composition similar to the cold-soil (warm-marginal) sites of our study area, where winters are still cold enough for frost but too warm for a persistent insulating snow cover.Furthermore, our results indicate that such reorganization will have implications for both the nitrogen cycling and soil microbial interactions, with subsequent consequences for plant growth.We suspect that the observed changes in our study were largely attributed to the negative effects of frost exposure rather than the positive effects of warming; hence, we cannot foresee strong changes in the microbiome due to warmer winters for sites with currently warm (non-frost affected) conditions.

F
I G U R E 3 Heatmap trees showing the taxonomic composition for Bacteria/Archaea (A) and fungi (C) as well as the core community of each (B, D).Only taxa with a summed (across all samples) proportion of reads >0.005 or 0.05 for bacteria/archaea (N = 90) and fungi (N = 85), respectively, are shown.In (A) and (B), tree nodes and branches are scaled and coloured to the relative abundance.In (B) and (D), only nodes and branches for members of the core community are coloured and labelled (N = 44 and 434 amplicon sequence variants [ASVs] for fungi and bacteria/archaea respectively).
Environmental characteristics.Note: Soil variables were measured in summer 2016 prior to the start of the experiment.Air temperatures (average, average minimum and absolute minimum temperature) are given as the average measurements from the experimental period: 1 December 2016 to 12 March 2017.'Ave.min.' and 'Abs.min.' refers to the statistical mean and the absolute minimum of the lowest air temperatures measured in the period.CN, the soils C:N ratio, OM, organic matter content (for measurement details on soil parameters see Site characteristics under Experimental Procedures).Sites are ordered according to soil temperature measurements during the experimental period from coldest to warmest from top to bottom (as shown in Figure1), which presents a different ordering than the air temperature variables in this table.The site abbreviation (and name) stand for the local forest municipality.
T A B L E 1 amplification products were purified using HighPrep PCR clean-up (MagBio Genomics, Gaithersburg, MD, USA) using a 0.65:1 (beads:PCR reaction) volumetric ratio.A second PCR reaction was performed to add Illumina sequencing adapters and sample-specific dual indexes (IDT Integrated DNA Technologies, Coralville, IA, USA) using PCRBIO HiFi (PCR Biosystems Ltd., London, UK) for 15 amplification cycles.Second PCR products were purified with HighPrep PCR clean-up as described for the first PCR.Sample concentrations were normalized using the SequalPrep Normalization Plate (96) Kit (Thermo Fisher Scientific, Waltham, MA, USA).The libraries were then pooled and concentrated using DNA Clean and Concentrator-5 Kit (Zymo Research, Irvine, CA, USA).The 16S and ITS2 pool concentrations were determined using the Quant-iT High-Sensitivity DNA Assay Kit (Life Technologies, Carlsbad, CA, USA).The libraries were denatured, diluted to 8 pM and sequenced following the manufacturer's instructions on an Illumina MiSeq platform using Reagent Kit v2 (2 Â 250 cycles; Illumina, San Diego, CA, USA) at the Section for Microbiology, University of Copenhagen.