Varying Vegetation Composition, Respiration and Photosynthesis Decrease Temporal Variability of the CO2 Sink in a Boreal Bog

We quantified the role of spatially varying vegetation composition in seasonal and interannual changes in a boreal bog’s CO2 uptake. We divided the spatially heterogeneous site into six microform classes based on plant species composition and measured their net ecosystem exchange (NEE) using chamber method over the growing seasons in 2012–2014. A nonlinear mixed-effects model was applied to assess how the contributions of microforms with different vegetation change temporally, and to upscale NEE to the ecosystem level to be compared with eddy covariance (EC) measurements. Both ecosystem respiration (R) and gross photosynthesis (PG) were the largest in high hummocks, 894–964 (R) and 969–1132 (PG) g CO2 m−2 growing season−1, and decreased toward the wetter microforms. NEE had a different spatial pattern than R and PG; the highest cumulative seasonal CO2 sink was found in lawns in all years (165–353 g CO2 m−2). Microforms with similar wetness but distinct vegetation had different NEE, highlighting the importance of vegetation composition in regulating CO2 sink. Chamber-based ecosystem-level NEE was smaller and varied less interannually than the EC-derived estimate, indicating a need for further research on the error sources of both methods. Lawns contributed more to ecosystem-level NEE (55–78%) than their areal cover within the site (21.5%). In spring and autumn, lawns had the highest NEE, whereas in midsummer differences among microforms were small. The contributions of all microforms to the ecosystem-level NEE varied seasonally and interannually, suggesting that spatially heterogeneous vegetation composition could make bog CO2 uptake temporally more stable.


ABSTRACT
We quantified the role of spatially varying vegetation composition in seasonal and interannual changes in a boreal bog's CO 2 uptake. We divided the spatially heterogeneous site into six microform classes based on plant species composition and measured their net ecosystem exchange (NEE) using chamber method over the growing seasons in 2012-2014. A nonlinear mixed-effects model was applied to assess how the contributions of microforms with different vegetation change temporally, and to upscale NEE to the ecosystem level to be compared with eddy covariance (EC) measurements. Both ecosystem respiration (R) and gross photosynthesis (P G ) were the largest in high hummocks, 894-964 (R) and 969-1132 (P G ) g CO 2 m -2 growing season -1 , and decreased toward the wetter microforms. NEE had a different spatial pattern than R and P G ; the highest cumulative seasonal CO 2 sink was found in lawns in all years (165-353 g CO 2 m -2 ). Microforms with similar wetness but distinct vegetation had different NEE, highlighting the importance of vegetation composition in regulating CO 2 sink. Chamberbased ecosystem-level NEE was smaller and varied less interannually than the EC-derived estimate, indicating a need for further research on the error sources of both methods. Lawns contributed more to ecosystem-level NEE (55-78%) than their areal cover within the site (21.5%). In spring and autumn, lawns had the highest NEE, whereas in midsummer differences among microforms were small. The contributions of all microforms to the ecosystem-level NEE varied seasonally and interannually, suggesting that spatially heterogeneous vegetation composition could make bog CO 2 uptake temporally more stable.

INTRODUCTION
Photosynthesis in boreal peatlands binds more CO 2 from the atmosphere than is released back through decomposition and plant respiration. For this reason, these ecosystems are sinks of atmospheric carbon (C) and especially important as a stable storage of 273-547 Pg of C (Turunen and others 2002;Yu 2011;Kleinen and others 2012), which equals to about one-third of soil C globally. However, in certain years the annual CO 2 balance may turn from uptake to release, because the C sink of boreal peatlands is small and sensitive to variations in temperature and moisture (Alm and others 1999;Waddington and Roulet 2000;Lund and others 2012) as well as light conditions (Nijp and others 2015). Due to the predicted increase in cloud cover and fluctuations in temperature and precipitation at northern latitudes (IPCC 2013), peatland C balance is likely to change. The fate of peatland C sink and storage under changing conditions has been simulated with process-based models (for example, Frolking and others 2010; St-Hilaire and others 2010; Wu and others 2011;Gong and others 2013). However, for such models more information is needed about the effect of spatially varying vegetation composition on the C cycling processes of peatland ecosystems. Within a peatland, distance of water table (WT) to the moss surface may vary widely, resulting in different microforms-hummocks, lawns and hollows-along the WT gradient, each being habitat for a particular plant community composition. Within-site WT variation is often especially pronounced in nutrient-poor boreal peatlands, that is, bogs, dominated by dwarf-shrub vegetation at the dry end of the WT gradient, gradually changing toward sedge-dominated, wet hollows. In the moss layer, Sphagnum species adapted to certain moisture conditions dominate at different points of the WT gradient (Hayward and Clymo 1982;Rydin 1993). These microforms differ greatly in net ecosystem CO 2 exchange (NEE), and one plant community may act as a C sink, whereas another is simulta-neously acting as a C source (Waddington and Roulet 2000).
Both ecosystem respiration (R) by plants and decomposers and gross photosynthesis (P G ) by plants have been found to be enhanced by the thicker aerobic peat layer in dry bog microforms, supporting more photosynthesizing vascular leaf area and aerobic decomposition. However, the spatial pattern of NEE is not straightforward; drier communities have been reported to have larger (Waddington and Roulet 2000), similar (Alm and others 1999;Bubier and others 2003a) and smaller (Strack and Waddington 2007) NEE as wetter ones. In previous studies, vegetation within the WT gradient is typically divided into two or three microform classes. Studying the CO 2 exchange using a finer classification could reveal important thresholds for R, P G and NEE and thus help explain the contrasting spatial patterns previously observed in bogs.
Peatland plant communities differ in their phenology, and consequently, in their seasonal development of R and P G . In spring, when the deciduous sedges have low leaf area, hummocks dominated by Sphagnum mosses and dwarf shrubs already are photosynthesizing (Moore and others 2006;Riutta and others 2007). Sedges in turn have their highest photosynthesis during their midsummer peak in leaf area (Korrensalo and others 2017). The seasonal changes in R are controlled by WT (Fenner and Freeman 2011) and temperature (Lafleur and others 2005), which differ in their importance in dry and wet microforms (Maanavilja and others 2011). In addition, when moisture and temperature regimes change between years, plant communities differ in their responses; for example, in a dry year, R, P G and NEE of some plant communities may decrease, while increasing or staying similar in other communities (Bubier and others 2003b). The seasonal timing and magnitude of changes in temperature and precipitation have been found to be an important control for the interannual variation of bog C uptake (Waddington and Roulet 2000;Lund and others 2012;Peichl and others 2018). Studying the C sink of boreal bog vegetation communities over several years could reveal which seasonal patterns occur persistently and which occur as a result of interannual changes in environmental conditions, for example, only during warm or cold years. However, to our knowledge, there are only two previous studies investigating the spatial variation in bog CO 2 exchange over more than two growing seasons (Strack and Waddington 2007;Pelletier and others 2011). Neither of these studies covers the whole growing season despite Varying Vegetation Composition, Respiration and Photosynthesis Decrease Temporal Variability the potential importance of spring and autumn photosynthesis for bog CO 2 uptake (Moore and others 2006;Korrensalo and others 2017).
In this study, we quantify the roles of plant communities with strongly differing species composition and WT depth for the CO 2 sink of a boreal bog over three growing seasons. To fulfill our aim, we (1) quantify the response of community and ecosystem-scale NEE to variations in light, temperature, moisture and leaf area, (2) determine how plant communities differ in their seasonal timing and magnitudes of P G , R and NEE and (3) examine variations in the ecosystem-level CO 2 sink by upscaling community-level NEE to an ecosystem-level estimate to be compared with ecosystemscale CO 2 flux measured by the eddy covariance method.

Study Site
The study site is an ombrotrophic bog (61°50¢N, 24°12¢E), located within the Siikaneva peatland complex in Ruovesi, Western Finland, in the Southern Boreal vegetation zone (Ahti and others 1968). The annual temperature sum in the area is 1318 degree days, annual rainfall is 707 mm, onethird falling as snow, and the average annual, January and July temperatures are 4.2, -7.2 and 17.1°C, respectively (30 year averages from the Juupajoki-Hyytiä lä weather station, 10 km from the site).
The large spatial variation in the WT level at the site manifests as a mosaic of vegetation communities differing in their species composition (Table 1). Dwarf shrubs (mainly Calluna vulgaris) and few small Scots pines (Pinus sylvestris) dominate the vascular vegetation in the dry end of the WT gradient. Toward the deeper WT, the dominance of field layer vegetation is shifted from dwarf shrubs to ombrotrophic sedges (Rhynchospora alba, Carex limosa) and the herb Scheuchzeria palustris. Groundlayer vegetation varies from Sphagnum fuscum dominating the high hummocks to lawns covered mainly by S. rubellum and S. papillosum and to S. majus-dominated hollows. The site also has bare peat surfaces without mosses but a sparse cover of R. alba and open water pools without any vegetation. An eddy covariance (EC) tower was located in the center of the site.

Spatial Sampling
The variation in vegetation composition was disaggregated into seven microform classes (Table 1), open water and six vegetated microforms, mainly based on the areal cover of Sphagnum moss species, which are known to vary in relation to WT depth (Hayward and Clymo 1982;Rydin 1993). Eighteen permanent sample plots were established around the EC tower in three groups, one sample plot in each of the three groups. To examine how well the division into classes represents the variation in vegetation within the site, a systematic vegetation inventory was made within 30 m radius from the EC tower in July 2013. The inventory grid consisted of a total of 121 points (Appendix 1), where the percent cover of each species was visually estimated, and based on those covers, the microform class of each point was defined. The areal cover of microforms based on vegetation inventory is reported in Table 1. We used detrended correspondence analysis (DCA) to examine how plant species composition is structured and related to water table. Canonical correspondence analysis (CCA) was used to test how well the division into microform classes explains the variation in vegetation. The cover of each species was also estimated at the 18 permanent sample plots, and these data were included in both analyses to examine how well they represent the variation in vegetation within their microform class and within the whole study area. The analyses were made with Canoco 5 software (ter Braak and Š milauer 2012).

Leaf Area Index Measurements
To follow the development of photosynthesizing leaf area at the 18 permanent sample plots over the growing seasons, we measured seasonal development of vascular green area (VGA) (Wilson and others 2007). The number of leaves of each vascular plant species in each sample plot was counted biweekly between May and September in 2012-2015. To define the VGA of each species inside the plots, the number of leaves was multiplied by the area (m 2 ) of an average leaf on each measurement day. Average leaf size for each species was defined on each measurement day with a scanner from a sample of 10-15 leaves collected next to the plots. For Calluna vulgaris, Betula nana, Empetrum nigrum and Vaccinium oxycoccos, we measured the total length (cm) of stem containing green leaves within the plots instead of the number of leaves and average leaf area (m 2 ) per stem centimeter instead of average leaf size.

Net Ecosystem Exchange Measurements
We used the closed chamber method (for example, Riutta and others 2007) to measure the exchange of CO 2 between the vegetated land surface and the atmosphere in different environmental conditions. The measurements were conducted weekly or biweekly during the snow-free period in 2012-2014 and every three weeks in 2015. In each measure-ment, transparent plexiglass chamber (56 9 56 9 30 cm) was sealed airtight by placing it into the water-filled groove of an aluminum collar permanently installed in the sample plots. The CO 2 concentration, photosynthetic photon flux density (PPFD) and the temperature inside the chamber were recorded every 15 s using radiation and temperature sensors and an infrared gas analyzer (EGM-4, PP systems, UK). A cooling system and a battery-operated fan maintained the chamber temperature within 2°C of the ambient temperature during the measurement. At every plot, 3-4 measurements of 90-180 s were conducted each day: first in full light, then one or two times with PPFD level reduced by 40-90% with a mosquito net, and finally with an opaque shroud. WT level was measured in perforated plastic tubes installed next to each plot, and the temperature was measured 5 and 15 cm below the moss surface. Both WT and temperature were recorded simultaneously with the other measurements at each plot. NEE was calculated from the linear change in CO 2 concentration in the chamber headspace and expressed as mg CO 2 m -2 h -1 . CO 2 exchange measured in the dark represented the sum of autotrophic and heterotrophic respiration, R. Here, we use the sign convention where both R and P G get positive values, and positive NEE (P G minus R) indicates a sink of CO 2 from the atmosphere.

Response Modeling
To reconstruct the NEE flux of the six measured microforms and to compare their responses to environmental variables, we used a nonlinear mixed-effects model with the hyperbolic light saturation curve (Larcher 2003): where NEE trsij is the observed net CO 2 exchange per square meter and the predictor PPFD trsij is the photosynthetic photon flux density (lmol m -2 s -1 ) for measurement j on plot i on plot group s on Julian day r in year t. The parameters to be estimated are respiration (R trsi ), maximum rate of light-saturated photosynthesis (Pmax trsi ) and the level of PPFD where half of Pmax was reached (k trsi ). The residual error (e trsi ) is normally distributed with mean zero and constant variance. In the full model, R trsi , Pmax trsi and k trsi were written as linear functions of fixed predictors and random effects. Possible fixed predictors were VGA, air temperature, WT depth and microform. To account for shading of vascular plants at high VGA Varying Vegetation Composition, Respiration and Photosynthesis Decrease Temporal Variability values (Riutta and others 2007), VGA was transformed as VGA B = 1exp(-VGA). Air temperature was transformed as (max(T, v)) 2 , where a second-order polynomial response was assumed for temperatures above v and no effect for temperatures below v.
The final submodels for the parameters of NEE were ðmaxðT; vÞÞ 2 ; n t ; n tr ; n trs ; n trsi Þ ð3Þ where C krsi is the categorical predictor of microform (6 levels). T is the air temperature (°C) using a transformation v = 17 for R trsi and v = 12 for Pmax trsi . The four last terms represent the random effects. The model was fitted using the function nlme of the R program package nlme (Pinheiro and Bates 2000). A more detailed description of the response modeling is presented in Appendix 2.
Using the fixed part of the model [equations (1)-(4)], we reconstructed daily (g CO 2 m -2 d -1 ) and seasonal (g CO 2 m -2 growing season -1 ) cumulative P G , R and NEE fluxes at the 18 measurement plots during growing seasons 2012-2014 between May and September (Julian days 121-273), hereafter referred to as growing season. For the reconstruction, we used half-hourly PPFD and T data from SMEAR II measurement station 10 km from the site, and daily VGA of each plot calculated using the log-linear VGA development model (Wilson and others 2007). The model parameter estimates, including those of the random part, are presented in Appendix 3.

Model Validation
To validate the fixed part of our nonlinear mixedeffects model, we used the NEE measurements of growing season 2015 as an external validation dataset, which was not used for model construction. The correlation between the measured and reconstructed NEE in 2015 was explored using equations (1-4) without the random effects.

Upscaling and Comparison of Microforms
To upscale the modeled fluxes to an ecosystem-level NEE, P G and R estimates, we calculated the daily mean cumulative flux (g CO 2 m -2 d -1 ) for each microform during the growing seasons 2012-2014 as an average of the predicted daily cumulative fluxes in the three replicate plots. Ecosystemlevel daily flux was obtained as a weighted average of these microform-scale fluxes using two estimates of the contributions of microforms to the cumulative EC tower footprint as weights (Table 2). These contributions were obtained from two footprint models following Kormann and Meixner (2001) and Kljun and others (2015). Seasonal cumulative NEE, P G and R are reported in Table 2 as a mean of these two upscalings as well as 95% confidence intervals of these values based on the standard errors of the fixed part parameters (equations 1-4, Appendix 3). Open water without vegetation was found to have R close to zero (27.2 ± 23.5 mg CO 2 m -2 h -1 , mean ± standard deviation, n = 37) in eight small measurement campaigns conducted between May and September 2013, and thus the areal cover of water surface was included in the upscaling without a flux value by assuming that it had R and P G of zero. The small area covered by boardwalk was treated in a similar way in the upscaling. The upscalings were done for an area within about 60 m for the Kljun and others (2015) model and 130 m for the Kormann and Meixner (2001) model from the EC tower, estimated to be the radii within which 80% of the EC flux originated (see description below).
We compared the reconstructed seasonal R, P G and NEE (18 values of each process in each year) among years and microforms with a linear mixedeffects model with year and microforms as well as their interaction as potential fixed effects and plot group as a random effect. We used a conditional F test to test whether the fixed effects or their interaction significantly improved the model. For the fixed effects included in the model, we obtained pairwise comparisons by fitting the full model multiple times with a different year or microform acting as a default level, to which the other years and microforms were compared. This analysis was performed with function lme of the R package nlme (Pinheiro and Bates 2000).

Eddy Covariance Measurements
The upscaled daily and cumulative growing season NEE fluxes measured by the closed chamber method were compared with the NEE estimates obtained by the EC technique. The EC method offers an independent and direct estimate of the ecosystem-level exchange of energy and matter by measurement of vertical turbulent fluxes (for  (Mammarella and others 2016). The processing of raw data (including despiking, coordinate rotation, de-trending, time lag adjustment, spectral correction, density correction and spectroscopic correction for LI-7700) was performed following the widely accepted routines (Aubinet and others 2012). The fluxes were quality-controlled and filtered by sensor signal strength and low nighttime turbulence conditions using a threshold value of 0.1 m/s for the friction velocity (u * ). According to footprint calculations with Kormann and Meixner (2001) and Kljun and others (2015) models, over 80% of the EC source area is contained within a 130 m (Kormann and Meixner 2001) and 60 m (Kljun and others 2015) radius around the EC tower in most conditions. The NEE time series for the growing seasons 2012-2014 was gap-filled with a model (NEE = GPP -Reco). First, a respiration model was derived from the nocturnal measurements in the form Next, the GPP was calculated as NEE-R e and modeled with a function where the drivers include PPFD (taken from the SMEAR II station) and the annual footprint-scale VGA, with b = const describing the moss leaf area. As in the response models of chamber data, a transformation (1exp(-VGA)) was used to model self-shading. The resulting GPP is expressed in mg (CO 2 ) m -2 h -1 .

Site Conditions
The most important gradient in species composition across the site (DCA axis 1, eigenvalue = 0.647) was related to moisture, along which the order of the plots was relative to their WT ( Figure 1A, B). Vegetation formed a continuum along this gradient from high hummocks (HHU) to wet hollow (HO) and bare peat (BP) surfaces ( Figure 1B). At the wet end of the moisture gradient, there was more variation also along a second gradient ( Figure 1B, DCA axis 2, eigenvalue = 0.209) that separated HO and lawn (L) surfaces with almost 100% Sphagnum cover from BP without mosses. The division of plots into microform classes significantly explained the variation in species composition (CCA; pseudo-F = 32.1, p = 0.002). In all years, WT separated the microforms into three groups, and maximum growing season VGA was closely related to that division. Those groups were (1) HHU with the lowest WT and highest VGA, (2) hummocks (HU) and high lawns (HL) with intermediate WT and VGA, and (3) L, HO and BP with the highest WT ( Figures 1A, 2). In the plant communities of the third group, dominating deciduous sedges made the seasonal VGA development sharper than at other communities, especially in hollows ( Figure 2).
Year 2012 was the coldest, cloudiest and wettest with a mean growing season air T of 13.0°C, daytime PPFD of 627 lmol m -2 s -1 (Figure 2) and WT 3 cm below the moss surface (weighted average of the microform WTs in Figure 1A). Growing season 2013 was the warmest, with a mean daily temperature of 14.5°C. Water table depth and the amount of light were similar in growing seasons 2013 and 2014 ( Figures 1A, 2), with WTs and daytime mean PPFDs of -6.7 cm and 654 lmol m -2 s -1 in 2013 and -6.5 cm and 655 lmol m -2 s -1 in 2014. In 2012-2013, daytime PPFD was the highest in the beginning and middle of the growing season, whereas in 2014, daytime PPFD showed more variation than in the previous years (Figure 2).
VGA was the highest during the warmest year (2013; Figure 2), with a site-level maximum VGA of 0.50 calculated as an average of the microform VGAs weighted by their cover within the site (Table 1). Despite the higher amount of light and deeper WT in 2014 than in 2012, VGA was rather similar in those two years (Figure 2). In 2012, the mean WT was above the surface at the three wettest microforms ( Figure 1A), with low VGA at lawns and bare peat surfaces, but not at hollows in that year.

CO 2 Exchange in Different Microforms
In the evaluation of the chamber model [equation (1-4)], the measured fluxes in the independent dataset of 2015 agreed well with the modeled fluxes (r 2 = 0.83, slope = 0.78, intercept = -5.5, Appendix 4). When using the model for flux reconstruction of years 2012-2014, significant dif-ferences were found among plant communities in reconstructed seasonal R (F = 1211.12, p < 0.001) and P G (F = 308.61, p < 0.001). The reconstructed R for the whole growing season decreased from the driest high hummocks to the wettest bare peat surfaces, and only lawns and hollows did not significantly differ in their seasonal R (Figure 3). There was also a larger division of plant communities into two groups, formed by the three driest and three wettest surfaces with higher and lower seasonal R, respectively (Figure 3). In reconstructed seasonal P G , the differences along the WT gradient were similar to R, except that bare peat surfaces had lower photosynthesis than the two other plant communities in the wet end of the gradient (Figure 3).
The relative differences in reconstructed seasonal NEE values among microforms were smaller than in reconstructed R and similar to reconstructed P G values (Figure 3). Significant differences were found among four groups of plant communities (F = 103.16 and p < 0.001), which were from the smallest to largest seasonal NEE: BP, HHU&HU, HL&HO and L (Figure 3). In all years, L communities were the largest CO 2 sinks, with mean annual NEE varying from 194 to 287 g CO 2 m -2 growing season -1 (Figure 3).

Varying Vegetation Composition, Respiration and Photosynthesis Decrease Temporal Variability
Interannual variation was higher in reconstructed P G than R (Figure 3). For all three variables-R, P G and NEE-the seasonal reconstructed fluxes were the lowest in the coldest, cloudiest and wettest year 2012 and highest in the warmest, sunniest and driest year 2013 (Figure 3). The interaction between years and plant communities was not significant in R (F = 0.58, p = 0.821), P G (F = 0.50, p = 0.877) or NEE (F = 0.49, p = 0.882).

Within-Season Variation in the Microform-Scale CO 2 Exchange
For most of the growing season, plant communities kept their ranking in terms of reconstructed daily R and P G (Figure 5). This was not the case with NEE, the net of R and P G , where the absolute differences in NEE were smaller and the ranking varied over the growing seasons ( Figure 5). In every spring and autumn, L communities had the highest NEE ( Figure 5). In midsummer, plant communities had very similar NEE, except in the coldest and wettest year of 2012 when HL and HO had the highest midsummer NEE ( Figure 5). In 2012, the different timing of maximum VGA between dry and wet communities caused the three driest plant communities to have also their seasonal NEE maximum earlier than wetter communities (Figures 1, 5). BP had much lower (mostly negative) NEE than the other communities ( Figure 5).
The measured NEE fluxes followed the seasonal patterns of VGA development and air temperature (Figures 1, 4), and accordingly, the final NEE response model [equations (1-4)] included the fixed predictors of VGA B and Tair, which predicted both parameters P G and R. Of these parameters, P G had a higher threshold value v above which the temperature response is polynomial, and thus the reconstructed P G was more sensitive to Tair than reconstructed R (Figures 1, 5). This was also the case with VGA B , which had a stronger effect on P G than on R in the model (Appendix 3). As a result, R varied less seasonally than P G ( Figure 5). In 2012 and 2013, the especially high midsummer VGA peak in hollows (Figure 1) resulted in their P G exceeding that of lawns in the middle of the summer ( Figure 5). As PPFD is the main model driver (equation 1), short-term seasonal fluctuations in daytime PPFD were also clearly seen in the reconstructed daily P G (Figures 1, 5).

Interannual Differences in the Ecosystem-Scale CO 2 Exchange
The ecosystem-level NEE for the whole growing season was the lowest, 71 g CO 2 m -2 , in the wettest and coldest year 2012 and more similar between 2013 (148 g CO 2 m -2 ) and 2014 (114 g CO 2 m -2 ) ( Table 2). Lawn communities, with the highest areal cover within the site at 21.5%, had a high contribution to the annual NEE in all years,  Table 1, and species' full names in Appendix 5. 55-78%. Bare peat communities had similar WT to hollows (Figure 2) but were net sources of CO 2 in all years (Table 2). Both ecosystem-level R and P G were the highest in the driest year 2013, but R varied much less among the years than P G ( Table 2).

Comparison of the Upscaled CO 2 Exchange with Eddy Covariance Measurements
When upscaled to the ecosystem level, the daily NEE fluxes were generally smaller than NEE measured by the EC tower ( Figure 6). The cumulative  Table 1. growing season chamber-based NEE fluxes were less than half of those measured by the EC tower in all years (Table 2). Cumulative growing season R and P G were also lower than the values derived from EC measurements (Table 2). However, the ranking in cumulative growing season NEE and P G was similar among the years in both methods so that year 2012 had the lowest and year 2013 the highest estimate (Table 2). According to both methods, there was much less interannual variability in ecosystem-scale R than in P G . The magnitude of interannual variation in NEE was smaller in upscaled chamber values (71-148 g CO 2 m -2 growing season -1 ) than in the values based on EC measurements (173-405 g CO 2 m -2 growing season -1 ) ( Table 2).

DISCUSSION
Over the growing season, all studied microforms except BP were net sinks of CO 2 from the atmosphere over the studied years, and L communities with an intermediate WT had the highest cumulative seasonal NEE, that is, CO 2 sink, at our site. Similarly, in a poor fen, lawns had the highest NEE (Strack and others 2006), whereas in two earlier studies on bogs (Alm and others 1999; Moore and others 2002) the spatial differences in measured NEE were very small. However, our results contradict several other studies in different types of peatlands, where hollows have been found to be smaller C sinks than drier microforms (Waddington and Roulet 2000;others 2006, 2007;Strack and others 2006;Riutta and others 2007). On the other hand, the spatial variation in R and P G at our site gains support from previous studies; hummocks had higher R (Alm and others 1999; others 2006, 2007;Strack and others 2006) due to a thicker aerobic peat layer, and higher P G (Laine and others 2007; Munir and others 2014) due to a higher VGA than wet microforms. Many earlier studied sites have been more homogeneous, without a thick dwarf-shrub cover at the dry end of the WT gradient (Riutta and others 2007), or bogs where the hollows dominated by sedges and hollow Sphagna are either lacking (Lafleur and others 2005;Bubier and others 2003a) or replaced by open water pools (Waddington and Roulet 2000). Our results highlight the importance of heterogeneous plant species composition in regulating the spatial variability of the CO 2 sink. With the fine-scale, vegetation-based microform classification used in this study, we observed two Varying Vegetation Composition, Respiration and Photosynthesis Decrease Temporal Variability pairs of microforms, HU&HL and L&HO, that had very similar WT but distinct vegetation and different cumulative CO 2 sinks in all years. In addition, in our response model microform class was a better predictor for NEE than temporally varying WT.
The high NEE of lawn communities resulted from their highest CO 2 sink of all microforms in spring and autumn, a seasonal pattern in NEE that persisted over the three measured years. Lawns do not have a high VGA in spring and autumn, but they do have a very high cover of Sphagnum species adapted to intermediate WT, known to be more efficient photosynthesizers than hummocks species (Gunnarsson 2005;Granath and others 2009;Laine and others 2011). Although Sphagna are less efficient photosynthesizers than vascular plants (Moore and others 2002;Leppä lä and others 2008;Otieno and others 2009;Korrensalo and others 2016), in ecosystems where VGA is low (< 1.4 m 2 m -2 ), as in our site, mosses may account for up to 96% of the CO 2 exchange (Douma and others 2007). Also, Sphagna have been found to have a strong role in the ecosystem-scale C sink in early spring, when vascular photosynthesis has not yet started (Moore and others 2006;Korrensalo and others 2017). Hollows, which have even more efficiently photosynthesizing Sphagnum species than lawns (Gunnarsson 2005;Granath and others 2009;Laine and others 2011;Korrensalo and others 2016), did not have the highest NEE in early and late summer. Hollows and lawns were rather similar in their WT, R and high Sphagnum cover. However, capitula of lawn Sphagna grow more densely than those of hollow species (Korrensalo and others 2018), which may cause the slightly higher spring and autumn P G of lawns. As observed before (Jä rveoja and others 2018), midsummer peak in NEE was strongly related to timing and magnitude of VGA development. Midsummer differences in NEE among plant communities were small, except for bare peat surfaces, which were mostly carbon sources.
Despite the low number of plant species, the functional diversity of bog ecosystems is high, meaning the presence of plant functional types (PFTs) with different physiology, morphology, seasonality and resource requirements (Tilman and others 1997;Cadotte and others 2008). Based on manipulation experiments (Ward and others 2009;Kuiper and others 2014) or observed natural variation (Alm and others 1999;Bubier and others 2003b), bog plant communities and PFTs are known to have differential responses to changes in environmental conditions. The functional diversity of bogs has been found to make their C sink more stable over a growing season than the C sink of fens, which often have functionally more homogeneous, sedge-dominated vegetation (Bubier and others 1998;Leppä lä and others 2008). Different responses of microforms with differing vegetation were also observed in this study, where contributions of plant communities to the ecosystem-scale CO 2 sink varied both within the growing seasons and between years with different temperature, WT and PPFD. Bog PFTs are known to respond differently to changing environmental conditions. Sphagna growing in the wet microforms respond with changed growth rate both to WT and to temperature variations more readily than hummock species (Robroek and others 2007). The abundance and biomass production of vascular PFTs are also known to react in opposite directions as a result of WT drawdown (Mä kiranta and others 2018), as deciduous vascular plants have generally less transpiration-decreasing morphological adaptations than evergreen dwarf shrubs. Over wide environmental gradients, bog species composition may change completely, but at the same time, the functional identity of bog ecosystems may remain rather similar (Robroek and others 2017). Nonetheless, more research is needed to determine whether the large diversity of PFTs and the varying responses of vegetation communities could, over timescales from several years to decades, make the C sink of bogs more stable than of peatlands with more functionally homogeneous vegetation composition, and which mechanisms would be critical for this stability. Studies reporting C fluxes for multiple years in varying types of peatlands (Peichl and others 2014; Helfter and others 2015; Lund and others 2015) could enable a multiyear comparison of peatland ecosystems of similar climatic conditions but different vegetation structure. Further, the variability of plant species and functional traits in bog ecosystems should be studied in relation to their C sink properties.
The site-level NEE measured with either EC or chamber methods fell within the range of several previously measured boreal bogs (Laine and others 2007;Koehler and others 2011;Petrescu and others 2015;Wilson and others 2016). It was lower than measured in a temperate bog with more shrubs and higher levels of PPFD others 2001, 2003), and higher than the NEE of 67.8-72.2 g CO 2 m -2 a -1 measured in a coastal blanket bog (Lund and others 2015). Because our site was rather wet for a bog, a deeper summertime WT increased the growing season CO 2 sink measured by EC tower from 173 to 405 g CO 2 m -2 growing season -1 between years 2012 and 2013, instead of decreasing it, as in drier bogs (Roulet and others 2007;Helfter and others 2015).
The different magnitude of fluxes estimated by chambers and EC tower raises questions about the possible error sources of the two methods. A mismatch between upscaled chamber fluxes and EC tower flux, although smaller than ours, was also found by Maanavilja and others (2011), where the high NEE peaks observed by the EC tower were not seen in the upscaled chamber fluxes. In chamber data, the choices made in modeling may cause different outcomes (Laine and others 2009). In our NEE model, WT was only included as a predictor in the form of microform, although WT is known to be an essential environmental factor explaining R, P G and NEE in peatlands (for example, Tuittila and others 2004;Laine and others 2006;Riutta and others 2007). However, the seasonal rhythm of WT variation did not differ among plant communities (p > 0.05 in all years, data not shown), and microform also includes the variation related to species composition. Additionally, our NEE model is supported by the model validation. In EC measurements, comparison of several towers at the same site has revealed that spatial variability causes substantial differences in the measured CH 4 flux magnitude and diel course (Peltola and others 2015), but as far as we know such analysis has not been conducted on CO 2 fluxes. Gap filling is an important source of uncertainty in cumulative EC fluxes (Moffat and others 2007), and most of this uncertainty originates from the different treatments of long gaps, which in our case occur mainly at the beginning and end of the growing season. Based on the uncertainty of the EC gap filling method used in this study, reported in Moffat and others (2007), the 95% confidence intervals for EC-based seasonal NEE values can be estimated to be approximately 137-209, 348-462 and 318-404 for growing seasons 2012, 2013 and 2014, respectively. In addition, measured EC fluxes are prone to systematic and random uncertainties, usually ranging between 10 and 30% of the cumulative seasonal flux (for example , Baldocchi 2003;others 2006, 2016). Considering the confidence intervals of EC-based and chamber-based NEE (Table 2), the uncertainty ranges of the seasonal values based on the two methods did not overlap. The difference between EC-based NEE and chamber-based NEE was the largest in 2013, when it was higher than the mean annual CO 2 sink reported for northern (arctic to temperate) open peatlands (Petrescu and others 2015). Further, in our data interannual variability in NEE was larger for EC data (173-405 g CO 2 m -2 growing season -1 ) than for upscaled chamber fluxes (77-141 g CO 2 m -2 growing season -1 ). The difference between the two methods was thus substantial for interpreting both the magnitude and the interannual patterns of NEE. Based on these results, the observed difference between the most widely used methods for NEE estimation warrants more research to reliably predict C sink of northern peatlands. This research should include more work on revealing the magnitude of error sources of the CO 2 flux measurements and data processing by both EC and chamber methods. In this process, these two methods will provide an invaluable comparison against each other's.

CONCLUSIONS
Microforms along the WT gradient of an ombrotrophic boreal bog had strongly varying species composition and large variation in their seasonal reconstructed P G and R. Species composition has the potential to regulate CO 2 sink strength of the microforms. The results show for the first time that two microforms with similar WT but distinct vegetation had a different seasonal NEE for three consecutive years. The seasonal reconstructed NEE was the highest in the most abundant microform, lawn. The seasonally and interannually varying contributions to ecosystem-scale NEE of microforms dominated by different plant functional types suggests this variability in vegetation composition may add to the temporal stability of the site CO 2 sink.
The reconstructed daily chamber fluxes were less than half of the EC fluxes, and the interannual variability in seasonal ecosystem-level NEE differed between the two methods. This result calls for exploring more the uncertainties related to both of the methods.

AC KNOWLEDGEMENT S
Open access funding provided by University of Eastern Finland (UEF) including Kuopio University Hospital. This work is supported by Faculty of Science and Forestry, University of Eastern Finland, Finnish Cultural Foundation and Academy of Finland (project code 287039), Academy of Finland Centre of Excellence (118780), Academy Professor projects (312571 and 282842) and ICOS Finland (281255 and 3119871). We are grateful for the help and facilities from Hyytiä lä Forest Research Station and thank Janne Sormunen, María Gutierrez, Franziska Rossocha and Laura Kettunen for the help in the field. We also thank Steve Frolking for revising the English language of the text and A. Lauré n for advices and opinions.

Co mpliance with E thical St andards
Conflict of interest The authors declare that they have no conflict of interest.
Data Avai labi lity The chamber and leaf area data associated with the paper are published in the PANGAEA repository under: https://doi.org/10. 1594/PANGAEA.905575 (Korrensalo and others 2019). Eddy covariance and environmental data are available in: https://avaa.tdata.fi/web/smart/s mear/download.

OPEN ACCESS
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4 .0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.