Wintertime Greenhouse Gas Fluxes in Hemiboreal Drained Peatlands

: The aim of this study is to estimate wintertime emissions of greenhouse gases CO 2 , N 2 O and CH 4 in two abandoned peat extraction areas (APEA), Ess-soo and Laiuse, and in two Oxalis site-type drained peatland forests (DPF) on nitrogen-rich sapric histosol, a Norway spruce and a Downy birch forest, located in eastern Estonia. According to the long-term study using a closed chamber method, the APEAs emitted less CO 2 and N 2 O, and more CH 4 than the DPFs. Across the study sites, CO 2 ﬂux correlated positively with soil, ground and air temperatures. Continuous snow depth > 5 cm did not inﬂuence CO 2 , but at no snow or a thin snow layer the ﬂuxes varied on a large scale (from − 1.1 to 106 mg C m − 2 h − 1 ). In all sites, the highest N 2 O ﬂuxes were observed at a water table depth of − 30 to − 40 cm. CH 4 was consumed in the DPFs and was always emitted from the APEAs, whereas the highest ﬂux appeared at a water table > 20 cm above the surface. Considering the global warming potential (GWP) of the greenhouse gas emissions from the DPFs in the wintertime, the ﬂux of N 2 O was the main component of warming, showing 3–6 times higher radiative forcing values than that of CO 2 ﬂux, while the role of CH 4 was unimportant. In the APEAs, CO 2 and CH 4 made up almost equal parts, whereas the impact of N 2 O on GWP was minor.


Introduction
Carbon dioxide (CO 2 ), methane (CH 4 ) and nitrous oxide (N 2 O) are three major greenhouse gases (GHG), which influence Earth's radiation balance. Since the industrial period, the amount of greenhouse gases in the atmosphere has increased due to human activities. According to WMO (World Meteorological Organization) [1], atmospheric greenhouse gas concentrations increased from 278 ppm in the pre-industrial period to 407.8 ppm for CO 2 , 722 to 1869 ppb for CH 4 and 270 to 331.1 for N 2 O in 2018, respectively.

Study Sites
The study was conducted at two abandoned peat extraction areas (APEAs) and two drained peatland forests (DPFs) in eastern and southeastern Estonia. The locations of the study sites are shown in Figure 1.
Atmosphere 2020, 11, x FOR PEER REVIEW 3 of 23 3 decrease as compared to the soil with a natural snow cover. The snow cover limits diffusion of CH4 from the atmosphere into the soil [53] and, therefore, thinner snow cover may enhance soil CH4 uptake.

Study Sites
The study was conducted at two abandoned peat extraction areas (APEAs) and two drained peatland forests (DPFs) in eastern and southeastern Estonia. The locations of the study sites are shown in Figure 1.  Supplementary Table S1). Peat depth was more than 50 cm. In all 4 study plots, drainage work had been carried out in the early 1970s. From October 2017 until April 2019, we measured greenhouse gas fluxes in two APEAs: Ess-soo (ES) and Laiuse (LA). At Ess-soo (57°54′53″ N, 26°41′51″ E) the replicate plots were on adjacent peat fields which were separated by ditches. The peat fields had a slightly convex surface, without trees. Hare's-tail cottongrass (Eriophorum vaginatum L.) coverage was 10 to 20%. At Laiuse (58°47′26″ N, 26°31′45″ E), the two plots were on adjacent peat fields which were separated by ditches. The most dominant plant species were E. vaginatum and bog haircap moss (Polytrichum strictum Bridel, J. Bot. (Schrader)). Pine and birch trees grew sparsely at the study site.
Both the Ess-soo and Laiuse sites are representative of Estonia's APEAs, and the Norway spruce and Downy birch sites are typical drained forests in Estonia's DPFs.  Supplementary Table S1). Peat depth was more than 50 cm. In all 4 study plots, drainage work had been carried out in the early 1970s. From October 2017 until April 2019, we measured greenhouse gas fluxes in two APEAs: Ess-soo (ES) and Laiuse (LA). At Ess-soo (57 • 54 53" N, 26 • 41 51" E) the replicate plots were on adjacent peat fields which were separated by ditches. The peat fields had a slightly convex surface, without trees. Hare's-tail cottongrass (Eriophorum vaginatum L.) coverage was 10 to 20%. At Laiuse (58 • 47 26" N, 26 • 31 45" E), the two plots were on adjacent peat fields which were separated by ditches. The most dominant plant species were E. vaginatum and bog haircap moss (Polytrichum strictum Bridel, J. Bot. (Schrader)). Pine and birch trees grew sparsely at the study site.
Both the Ess-soo and Laiuse sites are representative of Estonia's APEAs, and the Norway spruce and Downy birch sites are typical drained forests in Estonia's DPFs.

Meteorological Data
The meteorological analyses were based on air temperature measured hourly and calculated daily, as well as daily, monthly and annual precipitation from the Tartu Meteorological Station (58 • Figure S1).

Field Measurements
The static closed-chamber method [54] was used for the measurement of CO 2 , CH 4 and N 2 O fluxes. Closed chambers (conical, made of polyvinyl chloride (PVC), height-50 cm, Ø-50 cm, volume-65 L) were installed in four replicates per plot. The chambers were sealed with a water filled ring on the soil surface. In case there was snow cover, the chambers were placed on the snow. Snow was not removed inside the collars and the chambers were sealed as previously described. Gas sampling was carried out twice a month from October to April. Measurements consisted of five gas samples which were collected into previously evacuated 50 mL gas bottles for 1 h (at 0, 15, 30, 45 and 60 min). At each plot during each gas sampling session, the depth of groundwater table (cm) in the observation wells (Ø-50 mm, up to 1.5 m deep PVC) was determined, and the soil temperature was measured at four depths (10, 20, 30 and 40 cm) by a handheld Comet S0141 temperature logger with Pt1000TG8 sensors (Comet Systems Ltd., Rožnov pod Radhoštšem, Czech Republic). During each gas sampling session at each plot, groundwater parameters were measured from piezometers (pH, oxygen concentration, redox potential, temperature, electrical conductivity) by a handheld YSI Professional Plus Multiparameter Water Quality Instrument with a Quatro field cable (YSI Incorporated, Yellow Springs, OH, USA), and air and ground surface temperature by a handheld Comet S0141 temperature logger with Pt1000TG8 sensors (Comet Systems Ltd., Rožnov pod Radhoštšem, Czech Republic). In addition, at drained forests sites, the soil temperature (at 5 cm depth) was measured with temperature probes (model CS 107, Campbell Scientific Inc., Logan, UT, USA) and soil volumetric water content (VWC; depth-10 cm) was recorded with water content reflectometers (model CS615, Campbell Scientific Inc., Logan, UT, USA) at each plot. The automated abiotic data were stored as 1-h averages on a data logger (CR1000, Campbell Scientific Inc., Logan, UT, USA).

Soil Analyses
For soil analyses, 3 composite samples from LA, DB and NS and 8 composite samples from ES were analyzed, and each sample consisted of 24 sub-samples at the LA and ES sites and 3 sub-samples at the DB and NS sites (Table S2). The sub-samples were mixed and homogenized before laboratory analyses. In the homogenized soil samples, the content of dry matter, ash (3 h at 550 • C) and total C were estimated using a dry combustion method on a varioMAX CNS elemental analyzer (Elementar Analysensysteme GmbH, Langenselbold, Germany). The soil samples were also analyzed for total Atmosphere 2020, 11, 731 5 of 21 nitrogen according to Kjeldahl (Tecator ASN3313), NH 4 -N (Tecator ASN 65-32/84) and NO 3 -N (Tecator ASN 65-31/84). The ammonium lactate extractable solution, using the flame photometric method was applied to analyze the available K, Ca and Mg content (Tecator ASTN 9/84). The analyses were carried out at the Plant Biochemistry Laboratory of the Estonian University of Life Sciences.

Gas Analyses
The gas concentrations in the collected samples were determined using a Shimadzu GC-2014 gas chromatography system (Shimadzu GC-2014 with an electron capture detector, flame ionization detector, and Loftfield autosampler, Loftfield Analytics, Göttingen, Germany) at the laboratory of the Department of Geography, Institute of Ecology and Earth Sciences, University of Tartu, Estonia [55]. Emission rates for one plot were calculated as the average of five subsamples. Gas fluxes were calculated from the linear increase or decrease in gas concentration in each chamber with time, using a linear regression equation [56], and were corrected for air temperature according to the ideal gas equation. The following criteria were followed for quality control: R-squared (R 2 ) p-value < 0.05, R 2 > 0.77 (5 data points) or R 2 > 0.9 (4 data points) for CO 2 and R 2 p-value < 0.1, R 2 > 0.65 (5 data points) or R 2 > 0.8 (4 data points) for N 2 O and CH 4 . To achieve this, one data point was deleted of necessary. In case the difference between the minimum and maximum was smaller than 20 ppm/ppb (ppm applies to CO 2 and ppb to CH 4 and N 2 O), the R 2 value was not considered and the data were included in the analysis. In case of snow cover, the volume of chamber was recalculated considering snow depth and density.
For the comparison of the climate impact of GHGs, we calculated their global warming potential (GWP) according to the radiative forcing. In a 100-year timespan, the GWP for CH 4 and N 2 O is correspondingly 28 and 296 times larger than that of the equivalent amount of CO 2 [58].

Statistical Analyses
Principal component analysis (PCA) was performed on environmental parameters using the R package ade4 v. 1.7-15 [59]. Differences in PCA between the sites were evaluated using PERMANOVA with 9999 permutations. The pairwise comparisons were corrected with the Bonferroni method using the vegan v. 2.5-6 R package [60].
The normality of distributions was checked using the Shapiro-Wilk, Anderson-Darling, Kolmogorov-Smirnov, Lilliefors and Jarque-Bera tests. The distribution of gas data deviated from normal, and hence non-parametric tests were performed. The median, 25% and 75% percentile, and minimum and maximum values of variables are presented. We used the Kruskal-Wallis analysis of variance (ANOVA) test and Dunn's multiple comparison test to check the significance of differences between gas fluxes for different land-use categories and between different years at each study site, and the Spearman rank correlation to analyze the relationship between GHG fluxes and environmental parameters. p values were considered statistically significant after Benjamin-Hochberg correction. Statistical analysis was carried out using Statistica and XLSTAT. The level of significance of p < 0.05 was accepted in all cases, except in evaluating CH 4 and N 2 O flux regressions, when p < 0.1 was accepted.

PCA Analysis
According to the principal component analysis (PCA) of the soil, water and gas emission characteristics, the APEAs (ES and LA) were under similar conditions and significantly different from peatland forests (DB and NS) (p < 0.001 in all comparisons), which were in turn similar to each other ( Figure 2). The APEAs differed from the forests mainly in higher water table, dissolved O 2 Atmosphere 2020, 11, 731 6 of 21 and oxidation-reduction potential (ORP) values, and in lower pH and water electrical conductivity. Different temperatures were strongly negatively related to snow depth, but they showed no significant difference between the studied sites. Considering the gas fluxes, the peatland forests were characterized by higher CO 2 -C and N 2 O-N, and lower CH 4 -C emissions compared to the APEAs.
Atmosphere 2020, 11, x FOR PEER REVIEW 6 of 23 6 Different temperatures were strongly negatively related to snow depth, but they showed no significant difference between the studied sites. Considering the gas fluxes, the peatland forests were characterized by higher CO2-C and N2O-N, and lower CH4-C emissions compared to the APEAs. More detailed relationships between the GHG fluxes and environmental characteristics were detected in the correlation analysis. More detailed relationships between the GHG fluxes and environmental characteristics were detected in the correlation analysis.

Soil CO 2 Flux
The daily average soil flux of carbon dioxide varied between −1.1 and 106 mg CO 2 The highest values were measured in October and April and the lowest values were measured in February ( Figure 3). There was no significant difference between the four study sites (Kruskal-Wallis ANOVA test, multiple comparison of mean ranks). The median values of CO 2 -C fluxes were similar between the APEAs and DPFs, though the emissions were a little higher at the DPFs ( Figure 3). There was no significant difference in CO 2 fluxes between the years at any study site (Kruskal-Wallis ANOVA test) ( Figure 4).
Across all study sites, the emissions correlated positively with soil temperature at all four depths (10, 20, 30 and 40 cm) ( Figure 5) and air temperature ( Table 2). The emissions correlated negatively with snow depth ( Figure S2). In the DPFs, the correlations between CO 2 flux and soil temperature (at 5, 10, 20, 30, and 40 cm depth) were stronger, and especially high in the Norway spruce forest (Table 2). 10, 20, 30, and 40 cm depth) were stronger, and especially high in the Norway spruce forest ( Table 2).
In the APEAs, at LA site the fluxes of CO2 showed a significantly positive correlation with redox potential (ORP) ( Table 2). In the DPFs, the fluxes of CO2 were significantly and positively correlated with water temperature and pH (only in NS site; Table 2).
The fluxes of CO2 at the Downy birch DPF site showed a positive correlation with redox potential (ORP) and soil electrical conductivity.    In the APEAs, at LA site the fluxes of CO 2 showed a significantly positive correlation with redox potential (ORP) ( Table 2). In the DPFs, the fluxes of CO 2 were significantly and positively correlated with water temperature and pH (only in NS site; Table 2).
The fluxes of CO 2 at the Downy birch DPF site showed a positive correlation with redox potential (ORP) and soil electrical conductivity.
The average values of cumulative winter half-year soil flux of CO 2 -C at the APEAs (ES and LA sites) and drained forests (NS and DB) were 47.8, 41.9, 159.8, and 194.7 g C m −2 , respectively. The average winter-day fluxes were 199, 228, 761, 927 mg C m −2 d −1 . Wintertime CO 2 release from the APEAs (ES and LA sites) and drained forests (NS and DB) on average accounted for 18, 12, 21 and 20% of the total annual emission, respectively.

CH4 Fluxes
The average CH4-C emissions of varied between −97.7 and 1201.9 µg CH4-C m −2 h −1 . There was a significant difference between the APEAs and DPFs (Kruskal-Wallis ANOVA test, multiple comparison of mean ranks; p < 0.05), with the highest amounts emitted from the APEAs ( Figure 6). The negative values, which indicate methane consumption, were registered mainly from the DPFs. There was no clear temporal pattern, but in general, the fluxes were higher from October to January and lowest in February and March ( Figure 6). There was no significant difference in the CH4 fluxes between the different years at the APEAs and the Norway Spruce site (Kruskal-Wallis ANOVA test) (

CH 4 Fluxes
The average CH 4 -C emissions of varied between −97.7 and 1201.9 µg CH 4 -C m −2 h −1 . There was a significant difference between the APEAs and DPFs (Kruskal-Wallis ANOVA test, multiple comparison of mean ranks; p < 0.05), with the highest amounts emitted from the APEAs (Figure 6). The negative values, which indicate methane consumption, were registered mainly from the DPFs. There was no clear temporal pattern, but in general, the fluxes were higher from October to January and lowest in February and March ( Figure 6). There was no significant difference in the CH 4 fluxes between the different years at the APEAs and the Norway Spruce site (Kruskal-Wallis ANOVA test) (Figure 7). CH 4 fluxes in winters 2014/2015 and 2018/2019 showed statistically significant differences at the Downy birch DPF site (Kruskal-Wallis ANOVA test) (Figure 7).
At Laiuse APEA and the Norway spruce DPF site, the CH 4 fluxes correlated positively with water table depth ( Table 2). CH 4 flux showed the highest values at a water table >20 cm below the surface ( Figure S4). At the APEAs, the fluxes of CH 4 showed a significant (p < 0.05) positive correlation with soil temperature at all four depths (10, 20, 30 and 40 cm; Table 2). The fluxes of CH 4 showed a negative correlation with snow depth and a positive correlation with water table depth and water temperature at the Laiuse APEA site. The fluxes of CH 4 showed a positive correlation with ground surface temperature and air temperature at the Ess-soo APEA site. In the drained forests, the fluxes of CH 4 were significantly and negatively correlated with water temperature and soil temperature at all four depths (10, 20, 30 and 40 cm) ( Table 2). At the NS site, CH 4 correlated negatively with water electrical conductivity and positively with water table depth.
Average values of the cumulative winter half-year soil efflux of CH 4 -C in the abandoned peat extraction areas (the ES and LA sites) and drained forests (Spruce and Birch) were 0.25, 0.31, −0.11, and −0.07 g C m −2 , respectively. Average winter-day emissions were 1175, 1471, −177, and −140 µg CH 4 -C m −2 d −1 . Wintertime CH 4 release from the APEAs (ES and LA site) on average accounted for 31 and 52% of the total annual emission, respectively, and wintertime atmospheric CH 4 consumption in the drained forests (Spruce and Birch) on average accounted for 46 and (33%) of the total annual consumption, respectively. Even though in the birch forest methane consumption occurred in most winters, there was a slight methane release at the beginning and/or the end of winter (Figure 7b).
showed a negative correlation with snow depth and a positive correlation with water table depth and water temperature at the Laiuse APEA site. The fluxes of CH4 showed a positive correlation with ground surface temperature and air temperature at the Ess-soo APEA site. In the drained forests, the fluxes of CH4 were significantly and negatively correlated with water temperature and soil temperature at all four depths (10, 20, 30 and 40 cm) ( Table 2). At the NS site, CH4 correlated negatively with water electrical conductivity and positively with water table depth.

N2O Fluxes
The average emissions of N2O-N varied between −32.2 and 1440.8 µg N2O-N m −2 h −1 . There was a significant difference between the APEAs and DPFs (Kruskal-Wallis ANOVA test; multiple comparison of mean ranks; p < 0.05; Figure 8); the highest emissions in N2O fluxes between the years at the APEAs and Downy birch DPF site (Kruskal-Wallis ANOVA test; Figure 9). N2O fluxes in the 2014/2015 winter differed significantly from fluxes at the Norway spruce DPF site in the 2017/2018 and 2018/2019 winters (Kruskal-Wallis ANOVA test) (Figure 8). The highest peaks were registered in November and March, when freeze-thaw events were most common (Figure 9). Across all sites, emissions of N2O correlated negatively with water table depth. Thus, the highest N2O fluxes were observed at a water table depth of from −30 to −40 cm (Figure 10a). The N2O emission showed high values at soil temperatures around 0 °C and 8 °C (Figure 10b). N2O fluxes at the ES and LA sites were negatively correlated with air temperature. At the ES site, the N2O fluxes     (Table 2).

Differences in Peat Layer and Environmental Factors
There was a significant difference in peat parameters between the APEAs and DPFs (Figure 2). The top layer of peat had been removed from the peat extraction areas. The remaining peat layer consisted of compacted peat with likely recalcitrant organic carbon. Meanwhile, in the drained forests, only drainage had been conducted and no soil was removed. The most acidic soil occurred in the APEAs (especially the Ess-soo site), while the DPFs were the most neutral (pH values: 2.5, 3.0, 3.8, 4.0; in the Ess-soo and Laiuse APEAs and the Downy birch and Norway spruce DPFs, respectively). There were also significant differences in C/N ratio, which was considerably higher in the APEAs and lower in the DPFs: 45, 68, 15 and 19 in the Ess-soo, Laiuse, Downy birch and Norway spruce sites, respectively. Similarly, environmental factors such as precipitation, soil moisture, groundwater depth, air and soil temperature varied between the sites. These differences were crucial in the formation of the GHG fluxes in the sites (Appendix A).

Soil CO2 Fluxes
Average estimated wintertime CO2-C losses from the APEAs and drained forests were 42-48 and 169-195 g CO2-C m −2 , respectively. Wintertime CO2 efflux made up 10-25% of the annual release. This is similar to findings in earlier works, where wintertime CO2 fluxes have been estimated. The share of wintertime CO2 from annual respiration in a drained peatland forest was 21% [5], in a boreal agricultural peat soil 23% [36], in a minerotrophic fen 25% [61] and from a boreal forest 20% [31]. In general, winter CO2 emissions in temperate and boreal areas are in the range of 10 to 40% [61]. In Arctic areas with more than 200 days of snow cover, winter CO2 emissions are 10-30% of annual soil respiration [31]. The variation can be due to different soil types and ecosystems. In our study, the emission and percentage from annual release was higher in the DPFs and lower in the APEAs with sparse vegetation.
Our results closely correspond with those of other studies in which CO2 fluxes were well explained by air temperature and soil temperature [5,31,35,36,47,62]. Kim and Kodama [31] found strong correlation between wintertime CO2 flux and air temperature. Lohila et al. [36] have suggested that air temperature is the most important component that influences CO2 efflux throughout the whole year. Our results showed a slightly stronger correlation with upper layer soil temperatures than air temperature.
The correlations between CO2 emission and air and soil temperatures were stronger in drained forests. In all areas, the correlation was strongest with the topsoil (5-10 cm) temperatures and lower with deeper (40 cm) soil temperatures.
Waddington et al. [63] have pointed out the importance of the top layer (50 cm) of peat in CO2 production of peatlands, as this has higher substrate quality due to its proximity to organic matter inputs and better access to oxygen. This may also explain why in the APEAs (ES, LA) CO2-C In the spruce forest, the fluxes of N 2 O were positively correlated with water electrical conductivity and water table depth, and negatively correlated with soil temperature at 10 and 30 cm depth ( Table 2). The

Differences in Peat Layer and Environmental Factors
There was a significant difference in peat parameters between the APEAs and DPFs (Figure 2). The top layer of peat had been removed from the peat extraction areas. The remaining peat layer consisted of compacted peat with likely recalcitrant organic carbon. Meanwhile, in the drained forests, only drainage had been conducted and no soil was removed. The most acidic soil occurred in the APEAs (especially the Ess-soo site), while the DPFs were the most neutral (pH values: 2.5, 3.0, 3.8, 4.0; in the Ess-soo and Laiuse APEAs and the Downy birch and Norway spruce DPFs, respectively). There were also significant differences in C/N ratio, which was considerably higher in the APEAs and lower in the DPFs: 45, 68, 15 and 19 in the Ess-soo, Laiuse, Downy birch and Norway spruce sites, respectively. Similarly, environmental factors such as precipitation, soil moisture, groundwater depth, air and soil temperature varied between the sites. These differences were crucial in the formation of the GHG fluxes in the sites (Appendix A).

Soil CO 2 Fluxes
Average estimated wintertime CO 2 -C losses from the APEAs and drained forests were 42-48 and 169-195 g CO 2 -C m −2 , respectively. Wintertime CO 2 efflux made up 10-25% of the annual release. This is similar to findings in earlier works, where wintertime CO 2 fluxes have been estimated. The share of wintertime CO 2 from annual respiration in a drained peatland forest was 21% [5], in a boreal agricultural peat soil 23% [36], in a minerotrophic fen 25% [61] and from a boreal forest 20% [31]. In general, winter CO 2 emissions in temperate and boreal areas are in the range of 10 to 40% [61]. In Arctic areas with more than 200 days of snow cover, winter CO 2 emissions are 10-30% of annual soil respiration [31]. The variation can be due to different soil types and ecosystems. In our study, the emission and percentage from annual release was higher in the DPFs and lower in the APEAs with sparse vegetation.
Our results closely correspond with those of other studies in which CO 2 fluxes were well explained by air temperature and soil temperature [5,31,35,36,47,62]. Kim and Kodama [31] found strong correlation between wintertime CO 2 flux and air temperature. Lohila et al. [36] have suggested that air temperature is the most important component that influences CO 2 efflux throughout the whole year. Our results showed a slightly stronger correlation with upper layer soil temperatures than air temperature.
The correlations between CO 2 emission and air and soil temperatures were stronger in drained forests. In all areas, the correlation was strongest with the topsoil (5-10 cm) temperatures and lower with deeper (40 cm) soil temperatures.
Waddington et al. [63] have pointed out the importance of the top layer (50 cm) of peat in CO 2 production of peatlands, as this has higher substrate quality due to its proximity to organic matter inputs and better access to oxygen. This may also explain why in the APEAs (ES, LA) CO 2 -C emissions were lower in comparison with drained forests-Sphagnum moss had been removed from the sites and the area was covered with sparse vegetation. Though at the LA site there were more vegetation, the average water level was also higher (party flooded in spring during snow melt). Due to the high water level, the circumstances were unfavorable for further peat oxidation, which resulted in lower soil CO 2 efflux compared to the ES site. In both APEAs, the easier biodegradable upper layers of peat had been removed and the remaining deeper layers likely consisted of recalcitrant peat, which was not favorable for CO 2 production [64,65].
There was no clear relationship between CO 2 -C efflux and water level ( Figure S3). A significant difference in water table depth was found between the LA site and other sites (Kruskal-Wallis analysis of variance (ANOVA) test and Dunn's multiple comparison test, p < 0.01). The LA site was slightly flooded in springtime, up to 7 cm above the ground. CO 2 measurement in the presence of snow always provides challenges, especially in warming climate conditions [66]. Negative correlation with snow depth was found in all study sites (not statistically significant at ES site). The correlation with snow depth was stronger in the DPFs and especially high in one Downy birch forest microsite (R 2 = −0.75, p < 0.001). Soil temperature, which is affected by the thickness of snow and air-temperature, regulates the activity of microbes, which determine the amount of CO 2 flux from the soil [35]. If the ground is covered with snow, especially with deep snow, then the soils are wetter and warmer, which can increase heterotrophic respiration. In the case of less snow, the soils could freeze more deeply and result in lower CO 2 emissions. However, thin and sporadic snow cover may cause higher frequency of freeze-thaw events and winter CO 2 emissions could be larger with a warmer climate [67]. On the one hand, temporary ice cover formed on snow-free soil patches may store CO 2 in soil and degassing after soil is warming and ice is thawing [68], while on the other hand, freeze-thaw events may increase the soil-dissolved organic carbon (DOC) and dissolved organic nitrogen (DON) concentration [69] and, consequently, intensify the activity of soil decomposers [70] and GHG emission. The disturbance of fine roots can be a source for fresh available carbon for microorganisms [71,72]. Dense woody root mats in drained peatland forest sites (Table S1) are likely a reason why the CO 2 flux was higher here than that in the APEAs without any remarkable plant root system. Stielstra et al. [52] found that soil moisture is the primary determinant of carbon fluxes in seasonally snow-covered forest ecosystems. In this study, a weak positive correlation between CO 2 flux and soil moisture was found only in the spruce forest (NS).
At the APEAs, CO 2 consumption was registered in few measuring times. Generally, CO 2 flux is not negative out of the growing season. However, some studies showed that many evergreen plants and mosses maintain their photosynthetic activity throughout the winter [73][74][75]. This factor of CO 2 consumption can play a role in DPFs but not in APEAs with scarce vegetation.

CH 4 Fluxes
Average estimated wintertime CH 4 -C emissions from the APEAs and drained forests were 0.25-0.31 and −0.11-−0.07 g CH 4 -C m −2 , respectively. Wintertime CH 4 fluxes made up 31-52% of annual release at the APEAs and 33-49% of annual consumption in the drained forests. Average daily Atmosphere 2020, 11, 731 14 of 21 CH 4 fluxes from the APEAs were 1.2-1.5 mg CH 4 -C m −2 d −1 . Compared to previous studies in Minnesota peatlands, the emissions were lower (3-16 mg CH 4 m −2 d −1 ) [12]. However, wintertime fluxes made up 11-21% of annual emissions in the Minnesota peatlands, which is lower than our results. Our drained forests were annual sinks of CH 4 . This result corresponds well with earlier studies [76,77]. In the forest, the water table is relatively low due to canopy interception of precipitation and evapotranspiration [78], which may lead to a soil sink of CH 4 [76,77]. CH 4 consumption in our drained forests correlated negatively with soil temperature (strongest correlation with soil temperature at 20-40 cm) and positively with water level depth. Meanwhile, in our APEAs, the emissions of CH 4 correlated positively with air, water and soil temperature and water table depth. Low soil temperature slowed down CH 4 production in the APEAs, which is in accordance with earlier studies [5,13]. A part of methane emission in our APEAs was likely caused by the release of capped CH 4 from frozen surface layers: emission peaks occurred mostly when peat surface was slightly frozen. The fluxes were lower in mid-winter when snow cover was at its deepest. In Laiuse APEA, a significant correlation between CH 4 flux and snow depth was found. The N 2 O fluxes were influenced by soil properties, the ecosystem and the weather. No strong relationship was found between N 2 O emission rates and soil factors such as soil moisture (no correlation) and soil temperature (weak correlation). This result corresponds well with Flessa et al. [40]. The highest peaks of N 2 O emission occurred during freeze-thaw events. The months with the highest emissions were November and March. The high N 2 O peaks during the freeze-thaw events can be explained by a release of organic carbon as energy source for denitrification, by decomposing soil aggregates and killing soil organisms [10,33]. In the snowy winters, the N 2 O emissions were especially high in November and March. In mild winters with a thin snow cover or no snow at all, the emissions were higher also in mid-winter. In contrast, Brooks et al. [79] found that wintertime N 2 O fluxes were higher in the winters when snow cover formed earlier and lasted longer. Like Pärn et al. [80], who found an optimal range of soil water content for N 2 O emissions from organic soils, our study shows an intermediate water table depth (20-40 cm) at which the emission was highest (Figure 10a). The components of the determination of wintertime GWP were different for the APEAs and DBFs. In the DPFs, the flux of N 2 O was the main component, showing 3-6 times higher values than that of the CO 2 flux, while the role of CH 4 was of little importance. In the APEAs, CO 2 and CH 4 made up almost equal parts, whereas the impact of N 2 O on global warming was minor. In none of the study sites the slight consumption of CH 4 could compensate the GWP from CO 2 and N 2 O emissions.

Conclusions
Wintertime CO 2 and N 2 O fluxes in the DPFs were significantly higher than those from the APEAs. Vice versa, the mean CH 4 fluxes from the DPFs showed a cumulative consumption in the drained forests, whereas in the APEAs, emission dominated. Our hypotheses on snow cover effect was only partly supported, showing both a positive and negative impact on different sites and gases. In general, thin and scattered snow cover during the freeze-thaw periods initiates GHG fluxes. In our study, the winters were relatively mild and snow cover was rather thin. Across all study sites, CO 2 flux correlated positively with soil and ground and air temperature. A continuous snow depth of >5 cm did not influence CO 2 , while under no snow or thin snow fluxes showed the largest variation. Regarding the N 2 O and CH 4 fluxes, the correlation with snow depth varied between the gases and sites: in the APEAs and the spruce forest, snow cover slightly increased N 2 O flux, while CH 4 emission in the APEAs showed a negative correlation with snow depth. Most likely due to the freeze-thaw effect, the highest N 2 O emissions appeared at soil temperatures of around 0 • C, whereas the fluxes correlated negatively with soil temperature.
As a well-known phenomenon from previous studies, the highest N 2 O fluxes were observed at an intermediate water table depth of −30 to −40 cm, which corresponds with the optimal soil moisture level for microbial activity in the soil matrix.
Considering the global warming potential (GWP) of the greenhouse gas emissions from DPFs during the winters, the flux of N 2 O (GWP 265 times CO 2 equivalent) was the main component, showing 3-6 times higher values than that of the CO 2 flux, while the role of CH 4 (GWP 28 times CO 2 equivalent) was of little importance. In the APEAs, CO 2 and CH 4 made up almost equal parts, whereas the impact of N 2 O on global warming was minor.
In northern regions, milder winters are becoming more frequent due to the climate warming. An increase in freeze-thaw cycles will increase the likelihood of wintertime GHG emissions, especially from DPFs. Groundwater regulation around high levels in these areas would mitigate the emissions. Further analyses must focus on this and other possible measures to mitigate GHG fluxes from disturbed peatland systems. Our results might also be useful for the development of the GHG-based functional classification of soils [81].
Supplementary Materials: The following materials are available online at http://www.mdpi.com/2073-4433/11/7/ 731/s1, Table S1. Forest survey data in the Downy birch and Norway spruce sites. Average values are shown. Table S2. Soil physical-chemical parameters: Drained peatland forests-Downy birch (DB); Norway spruce (NS). Abandoned peat production areas: Ess-soo (ES), Laiuse (LA). Figure S1: Mean monthly precipitation (mm) (left) and monthly air temperature ( • C) (right) at Jõgeva (corresponding to Laiuse study site), Tartu-Tõravere (corresponding to the Spruce and Birch forest study sites) and Võru weather stations (corresponding to the Ess-soo study site) in the period of 2015-2019. Figure S2: The relationship between CO 2 -C efflux and snow depth (cm) at all study sites. Figure S3: The relationship between CO 2 -C efflux and water table depth (cm) at all study sites. Negative water table depth values show that the surface was flooded. Figure S4: The relationship between CH 4 -C emissions and water table depth (cm) at all study sites. Negative water table depth values show that the surface was flooded.