Temperature Control of Spring CO 2 Fluxes at a Coniferous Forest and a Peat Bog in Central Siberia

: Climate change impacts the characteristics of the vegetation carbon-uptake process in the northern Eurasian terrestrial ecosystem. However, the currently available direct CO 2 ﬂux measurement datasets, particularly for central Siberia, are insufﬁcient for understanding the current condition in the northern Eurasian carbon cycle. Here, we report daily and seasonal interannual variations in CO 2 ﬂuxes and associated abiotic factors measured using eddy covariance in a coniferous forest and a bog near Zotino, Krasnoyarsk Krai, Russia, for April to early June, 2013–2017. Despite the snow not being completely melted, both ecosystems became weak net CO 2 sinks if the air temperature was warm enough for photosynthesis. The forest became a net CO 2 sink 7–16 days earlier than the bog. After the surface soil temperature exceeded ~1 ◦ C, the ecosystems became persistent net CO 2 sinks. Net ecosystem productivity was highest in 2015 for both ecosystems because of the anomalously high air temperature in May compared with other years. Our ﬁndings demonstrate that long-term monitoring of ﬂux measurements at the site level, particularly during winter and its transition to spring, is essential for understanding the responses of the northern Eurasian ecosystem to spring warming.


Introduction
Boreal forests and peatlands are the major terrestrial biomes and large carbon (C) reservoirs in northern Eurasia, occupying 49% and 25% of Russia's land area, respectively [1,2]. Both ecosystems are considered essential C sinks in the global C cycle [2][3][4][5][6][7][8]. The role of annual or seasonal C-sink capacity in boreal forests can vary depending on temperature     1 29.7 10.1 Tower height (m) 2 29.4 9.8 Zero plan displacement 3 height (m) 13.4 1.675 for pine trees 0.335 for dwarf shrubs Roughness length (m) 4 3 0.375 for pine trees 0.075 for dwarf shrubs Leaf area index (m −2 m −2 ) 1-3.5 5 not available Soil type Podzol Histosol 1 Measurement height was taken as the height measured to the centre of the sonic anemometers. 2 Tower heights were taken as the distance from the ground fundament panels to the top of the tower plate. 3 Zero plan displacement height was computed as 0.67 × vegetation height. 4 Roughness length was computed as 0.15 × vegetation height. 5 Values obtained from the inventory and remote-sensing measurements [43,44].
The Zotino bog flux tower (hereafter ZB, 60 • 49 03 N, 89 • 23 20 E, 66 m a.s.l.) is situated approximately 2 km to the northeast of the ZOTTO. The average canopy height at the site was approximately 2.5 m, and the measurement height was 9.9 m (Table 1 and Figure 1). The type of peat is classified as an ombrotrophic [21]. The landscape was covered by a pine-dwarf, shrub-sphagnum (in Siberia, called ryam), hollow-ridge complex. They had a height of 0.4-0.6 m and were covered by plant communities consisting of dwarf pine (Pinus sylvestris f. litwinowii), which dominates the trees and dwarf shrubs (Chamaedaphne calyculata). The peat's calibrated age at the bottom of the bog ranged from 9397 ± 134 y BP at the edges to 13617 ± 190 y BP at the bog's center. The peat depth showed a wide range from 1.60 m to 5.10 m, increasing toward the bog's center.

Measurement System
Identical micrometeorological measurement systems were installed at the forest and bog sites ( Table 2). The EC system consisted of a three-dimensional ultrasonic anemometer (USA-1, Metek GmbH, Elmshorn, Germany) with integrated 55 W heating and closed-path infrared gas analyser (LI-7200, LiCor Biosciences, Lincoln, NE, USA) to measure CO 2 and H 2 O fluxes at 20 Hz frequency. An external diaphragm vacuum pump transported air to the gas analyser with a flow rate of 13 L min −1 at ambient atmospheric pressure. At the top of the towers, sensors measured the four radiation components and photosynthetically active radiation (PAR), air temperature (T a ), relative humidity, and atmospheric pressure. Soil or peat temperature (T s ) was measured by PT100 soil-temperature probes at six depths (0.02, 0.04, 0.08, 0.16, 0.32, and 0.64 m for forest and 0.04, 0.08, 0.16, 0.32, 0.64, and 1.28 m for bog, respectively). Soil moisture was measured by six sensors at both sites: two sensors at 0.08 m and one sensor at each depth of 0.16, 0.32, 0.64, and 1.28 m at the ZF and six sensors at a depth of 0.08 m at the ZB. Data collected from the EC system and meteorological measurements were stored on a data logger (CR3000, Campbell Scientific Inc., Logan, UT, USA). Details in setup of the EC system have been described by Park et al. [42]. Snow depths at both sites were measured manually at locations nearby the towers; however, the measurement intervals were irregular.

Post-Processing of Data and Quality Control
We applied the data processing scheme for the closed-path analysers mentioned in Mammarella et al. [45]. In the raw data post-processing step, spike detection, 2D coordinate rotation, and time lag adjustment as well as calculation of half-hourly turbulent fluxes were performed using the EddyUH software [45]. The thresholds of friction velocity (u*) for ZF and ZB were 0.2 m s −1 and 0.1 m s −1 , respectively. Details of the post-processing steps and quality control have been described by Park et al. [42]. The high-quality data available after the quality check and u* filtering led to a data coverage of 42-61% for ZF and 19-51% for ZB for 2013-2017. To calculate the cumulative CO 2 flux (NEE cum ), we gap-filled the half-hourly data using the REddyProc [46] in R [47]. Storage flux calculation was applied only for ZF data.

Data Anlaysis and Statistical Model
In this study, 'spring' is defined from the day of year (DOY) of 91-165 (01.04-14.06), 2013-2017. Daily CO 2 fluxes were summed for each 24-h period with gap-filled data. Regarding the EC convention, a negative NEE or CO 2 flux refers to net CO 2 uptake by the ecosystem, whereas a positive NEE indicates CO 2 release from the ecosystems to the atmosphere.
We used surface albedo (Alb) as a proxy for the status of snowmelt, as suggested in previous studies [22]. We determined the final day of snowmelt on which the 3-day moving average of Alb fell below 0.15 (15%) for the first time. To detect the final day of snowmelt, Shibistova et al. [28] used the daily mean soil temperature at a depth of 0.05 m and a 15-day moving average as the final day of snowmelt; however, in the present study, we used diurnal patterns of surface or peat temperature at a depth of 0.04 (T s04 ). To examine the effect of temperature accumulation on vegetation productivity, we defined the accumulated growing degree days for air temperature (CGDD Ta ) and surface soil or peat temperature (CGDD Ts04 ). We computed the cumulative sums of daily mean T a and T s04 above 0 • C (base temperature 0 • C) from DOY 99-165.
To identify the abiotic variables that are important in explaining the variability in daily CO 2 flux, we used the multivariate adaptive regression splines (MARS) regression model [48], specifically, the earth package [49] in R [47]. Daily CO 2 fluxes with corresponding abiotic variables were used as training data for both sites from DOY of 99 to 165 in 2013-2017 for ZF and the same days in 2014-2016 for ZB, respectively. Several previous studies have shown that spring CO 2 uptake is associated with four environmental factors: air temperature (T a ), T s04 , beginning day of snowmelt, and final day of snowmelt [21][22][23][24]. Therefore, we included T a , T s04 , and Alb to construct the MARS model. For ZF, we used nine abiotic variables as training data: PAR, Alb, T a , midday vapour pressure deficit (VPD), T s04 , and soil or peat temperature at depths of 0.08 m (T s08 ), 0.16 m (T s16 ), 0.32 m (T s32 ), and 0.64 m (T s64 ). For ZB, a total of seven abiotic variables were used as training data: the same as those for ZF but excluding T s08 and T s16 due to the long-term data gaps (>3 months) in the data for these two variables.
The MARS is a non-parametric regression method that can deal with both linear and nonlinear relationships and interactions between variables in the data using hinge functions [48]. Variable selection algorithms search the variables using both forward and backward stepwise selections. Variable importance is determined by a selection algorithm and is based on the number of model subsets (nsubsets), generalized-cross validation (GCV) score, and residual sum of squares (RSS). The GCV and RSS were scaled from 0 to 100. Higher GCV scores indicate variables with more explanatory power. The largest summed decrease of the RSS is scaled to 100, meaning that a large net decrease in the RSS is more important than a low RSS. Higher nsubsets means that variables are included in more subsets because they are important. Further statistical theory and application are described in detail at http://www.milbo.org/doc/earth-notes.pdf (accessed on 26 July 2021).
To investigate the impacts of cold spring weathers on the net ecosystem productivity (NEP), we used a rectangular hyperbolic light-response functional relationship [50] between PAR and NEP. We examined the differences in curve shapes and fit parameters for the year with the most frequent cold spring weather days (2014) and the year with the least frequent cold spring weather days (2015) with the following Equation (1): where PAR (µmol photon m −2 s −1 ) is the incident photosynthetically active radiation, A max (µmol CO 2 m −2 s −1 ) is the light-saturation point of CO 2 uptake, and α is the initial slope of the NEP-PAR (net CO 2 uptake at light saturation), and R d is the ecosystem respiration during the day (µmol CO 2 m −2 s −1 ). We used daytime data for which potential global radiation (R pot ) exceeded 20 W m −2 . To make it easier to visualize, we used NEP, which is opposite in sign to NEE. Three model parameters were estimated for the selected two periods estimated using the Levenberg-Marquardt method, implemented in the minpack.lm package [51] in R [47]. The Levenberg-Marquardt method is a widely used algorithm for analyzing non-linear light-response curves [52,53]. NEP sat was calculated with the obtained three parameters and by fixing PAR at 1500 µmol photon m −2 s −1 . The definition of NEP sat was the similar concept addressed by Migliavacca et al. [52] and Musavi et al. [54]. NEP sat helps to understand NEP at light saturation when PAR as 1500 µmol photon m −2 s −1 corresponds to light saturation.

Abiotic Controls of Spring CO 2 Fluxes
Temporal evolutions of daily CO 2 fluxes and associated abiotic drivers in 2016 showed that CO 2 flux variations generally corresponded to the rapid increasing T a and T s04 after snowmelt had completed ( Figure 2). We confirmed that these highly fluctuating features of fluxes and spring meteorological conditions were similar in other years (Supplementary Figure S1). CO 2 fluxes at ZF increased faster than at ZB after snowmelt. Similar features were observed in other years. Alb at ZB was remarkedly different from that at ZF because of the higher radiation absorption by forest. Sporadic warm spells temporarily led to a rapid net uptake of CO 2 at ZF. In contrast, net CO 2 uptake rates at ZB did not increase as rapidly as those at ZF because the Sphagnum peat was still under snow cover and therefore less affected by changes in T a . During the snowmelt period (DOY 127-145), the forest transitioned from a net CO 2 source to a net CO 2 sink. However, CO 2 fluxes were still highly variable during this time, fluctuating between net CO 2 source and net CO 2 sink. In a typical springtime, Ta was highly variable, ranging from −5 °C in April to 25 in May ( Figure 2). The amplitude of Ta was ~30 °C, which is typical of boreal spr weather. At the end of the snowmelt stage, the increase in Ta was considerable (by ~15 ° This may have led to the disappearance of snow cover on the ground surface, resulting the observed increase in Ts04 during this period to above 5 °C. As shown in Figure 2, also confirmed that Ts04 at ZB increased above 15 °C, which was ~5 °C warmer than at for every spring in the study period (Supplementary Figure S1).
At both ecosystems, the start of net CO2 uptake occurred during periods with cold frozen soil ( Figure 2). However, the magnitude of daily CO2 uptake was lower than − C m −2 d −1 until thawing of the surface soil (Ts04 < 0 °C) or snowmelt was complete. C uptake rates remained almost neutral and constituted a net CO2 sink only when Ts04 covered above 0 °C (DOY 120). As shown in Figure 2, data for both sites reveal that transition from net CO2 source to net CO2 sink can be disrupted temporarily by cold sp (DOY 110-135). This feature was more distinct in forest ecosystem. In a typical springtime, T a was highly variable, ranging from −5 • C in April to 25 • C in May ( Figure 2). The amplitude of T a was~30 • C, which is typical of boreal spring weather. At the end of the snowmelt stage, the increase in T a was considerable (by~15 • C). This may have led to the disappearance of snow cover on the ground surface, resulting in the observed increase in T s04 during this period to above 5 • C. As shown in Figure 2, we also confirmed that T s04 at ZB increased above 15 • C, which was~5 • C warmer than at ZF for every spring in the study period (Supplementary Figure S1).
At both ecosystems, the start of net CO 2 uptake occurred during periods with cold or frozen soil ( Figure 2). However, the magnitude of daily CO 2 uptake was lower than −1 g C m −2 d −1 until thawing of the surface soil (T s04 < 0 • C) or snowmelt was complete. CO 2 uptake rates remained almost neutral and constituted a net CO 2 sink only when T s04 recovered above 0 • C (DOY 120). As shown in Figure 2, data for both sites reveal that the transition from net CO 2 source to net CO 2 sink can be disrupted temporarily by cold spells (DOY 110-135). This feature was more distinct in forest ecosystem.
In 2016, the transition from net CO 2 source to net CO 2 sink occurred earlier in the forest than in the bog ( Figure 2). Although the soil was still frozen, if the temperature was favourable (~5 • C) for photosynthesis, the forest ecosystem transformed into a weak net CO 2 sink. For instance, the start of net CO 2 uptake at ZF was on DOY 114 in 2016. Mean T a from DOY 113 to 119 was 7.7 • C, which favoured the triggering of the start of net CO 2 uptake. Compared with ZF, ZB was a not yet an active CO 2 sink, probably because the ground layer peat mosses were partially covered by snow.
We used the MARS model to examine the relative importance of factors controlling CO 2 fluxes over the study period. The fitting of this model showed fairly good predictive performance (84-89%; Table 3). For ZF, the MARS model identified six variables (T s04 , PAR, Alb, T a , T s16 , and T s64 ) as the important influences on CO 2 fluxes out of the nine analysed. The percentage of variance in daily CO 2 fluxes explained by the model (R 2 ) was 84%. For ZB, the MARS model identified five variables (T s04 , Alb, T s64 , PAR, and T s32 ) as the important influences on CO 2 fluxes out of the seven variables analysed, with R 2 of 89% at ZB. Although we did not use ecological model, results show that the MARS model captured the controlling factors of CO 2 fluxes fairly well by considering non-linear relationships between modelled CO 2 fluxes and abiotic drivers.

Impacts of Cold Weather on Photosynthesis-Related Parameters
During the springtime, T a was highly variable, with several apparent cold and warm phases ( Figure 2). To examine the overall impacts of cold weather on CO 2 fluxes, we compared the photosynthesis-related parameters between the year having the most frequent cold spring weather days (2014) and the year having the least frequent cold spring weather days (2015) using the PAR-NEP relationship ( Figure 3 and Table 4). Here, we defined cold spring weather as the day with a minimum daily T a < 0 • C. For both sites, 2014 had the most frequent cold weather days (13 days for ZF and 15 days for ZB), and 2015 had the least frequent cold weather days (0 days for ZF and 2 days for ZB) compared with the entire study period (6 days for ZF and 10 days for ZB, 5-year mean values).  Table 4.
It is a typical feature that half-hourly CO2 fluxes during spring daytime are hig scattered due to large temperature changes or cloud cover ( Figure 3). Nevertheless, light-response curves of both ecosystems showed typical shapes of rectangular hyperb curve. The curve showed a rapid and linear increase in NEP at lower PAR levels, wi slow increase to reach maximum NEP at high PAR levels ( Figure 3). Bogs fixed CO much lower PAR levels compared with forests. Under the same environmental conditi for instance, needle-leaves at the top of the forest canopy are likely to warm sufficie to photosynthesise compared with the understory vegetation, such as that found at bog site. Due to the wider range of NEP for forest than those for bog, residual stand errors for ZF were larger than those for ZB (Table 4). Mean Ta values after snowme 2014 and 2015 were ~5 °C and ~13.5 °C, respectively. As a result, both ecosystems seem have higher Amax and NEPsat in 2015 than in 2014. Amax from the rectangular hyperb photosynthetic light-response function had a higher light saturation level than the from the non-rectangular hyperbola curve because the latter had a stronger flattenin the light-compensation point. Therefore, NEPsat better represents values of the maxim CO2 uptake rate of ecosystems. Compared with the bog, the forest had approxima three times higher Amax and NEPsat.
Quantum yields (α) at both ecosystems were similar in 2014 and increased by m than 50% in 2015 (Table 4). In this period, the light-response curves for ZB reached at li saturation point earlier than for ZF, indicating that the bog reached its maximum ab to fix CO2 earlier than for the forest due to lower light-use efficiency (Figure 3). Inter ingly, the relative changes in NEPsat of ZF and ZB in 2015 were approximately +32% 74%, respectively. Presumably, the contribution of soil microbes in the ecosystem resp tion term under the warm spring in 2015 may be greater than in 2014.  Table 4. Table 4. Model parameters of the rectangular light-response function for the year with the most frequent cold weather days (2014) and for the year with the least frequent cold weather days (2015) for both ZF and ZB. Numbers in parentheses denote the standard errors of the parameters. n is the number of half-hourly data. As NEP sat is a function of A max , α, R d , and PAR at saturation PAR of 1500 µmol photon m −2 s −1 , the standard deviations of NEP sat were taken as the averaged maximum and minimum standard errors of each parameter.

Site
Year It is a typical feature that half-hourly CO 2 fluxes during spring daytime are highly scattered due to large temperature changes or cloud cover ( Figure 3). Nevertheless, the light-response curves of both ecosystems showed typical shapes of rectangular hyperbola curve. The curve showed a rapid and linear increase in NEP at lower PAR levels, with a slow increase to reach maximum NEP at high PAR levels ( Figure 3). Bogs fixed CO 2 at much lower PAR levels compared with forests. Under the same environmental conditions, for instance, needle-leaves at the top of the forest canopy are likely to warm sufficiently to photosynthesise compared with the understory vegetation, such as that found at the bog site. Due to the wider range of NEP for forest than those for bog, residual standard errors for ZF were larger than those for ZB (Table 4). Mean T a values after snowmelt in 2014 and 2015 were~5 • C and~13.5 • C, respectively. As a result, both ecosystems seem to have higher A max and NEP sat in 2015 than in 2014. A max from the rectangular hyperbola photosynthetic light-response function had a higher light saturation level than the one from the non-rectangular hyperbola curve because the latter had a stronger flattening at the light-compensation point. Therefore, NEP sat better represents values of the maximum CO 2 uptake rate of ecosystems. Compared with the bog, the forest had approximately three times higher A max and NEP sat .
Quantum yields (α) at both ecosystems were similar in 2014 and increased by more than 50% in 2015 (Table 4). In this period, the light-response curves for ZB reached at light-saturation point earlier than for ZF, indicating that the bog reached its maximum ability to fix CO 2 earlier than for the forest due to lower light-use efficiency ( Figure 3). Interestingly, the relative changes in NEP sat of ZF and ZB in 2015 were approximately +32% and 74%, respectively. Presumably, the contribution of soil microbes in the ecosystem respiration term under the warm spring in 2015 may be greater than in 2014.

Interannual Variability of Spring Cumulative NEE cum
Overall variations in spring NEE cum were more distinct at ZF compared with ZB, as shown by daily CO 2 fluxes in Figure 2. Spring NEE cum between 2014 to 2016 ranged from −27.6 to −37.2 g C m −2 at ZF and from −7.7 to −14.9 g C m −2 at ZB ( Figure 4 and Table 5). Compared with the features at ZB, NEE cum at ZF showed larger amplitudes and higher uptake rates, which indicate that forest is a stronger spring CO 2 sink than the bog. When thick snow depth rapidly decreased with increasing temperature, between DOY 120 and 130, 2014-2016, both ecosystems released CO 2 in the atmosphere. Table 5. Cumulative net ecosystem exchange of CO 2 (NEE cum ; g C m −2 ), the day of year (DOY) of transition from NEE cum source to NEE cum sink (NEE cum,tran ; DOY), cumulative growing degree days of air temperature (CGDD Ta ; • C) corresponding to NEE cum,tran , cumulative growing degree days of surface soil temperature (CGDD T04 ; • C) corresponding to NEE cum,tran , and mean slope ( ∂NEE cum ∂DOY ; g C m −2 d −2 ) after reading peak of NEE cum for ZF and ZB for the period DOY 99-155 in 2014-2016. To compare the rate of change of interannual variability in NEE cum , the mean slopes of the curve ( ∂NEE cum ∂DOY ) following peak NEE cum were examined (Table 5). Both ecosystems showed steeper declining slopes after the peak of NEE cum in the three years studied, particularly the steepest in the warmest year (2015); ∂NEE cum ∂DOY was −1.62 g C m −2 d −2 for ZF and −0.92 g C m −2 d −2 for ZB. For both ecosystems, 2015 was the year with the largest NEE cum for both the 3-year period (2014-2016) and the entire study period. In 2015, T a in May was~10 • C, which was 4.2 • C higher than the 5-year mean (5.8 • C) ( Figure 5). For both ecosystems, NEE cum was strong and negatively correlated with mean T a in May (R 2 = 0.69, P = 0.081 for ZF and R 2 = 0.79, P = 0.110 for ZB). Unusually high mean T a in May corresponded to the highest NEE cum , meaning that the ecosystems were the strongest CO 2 sinks in the warmest year.  Rapid increasing trends in accumulated growing degree days of air and soil temperatures were followed by the steep decreasing slope in NEE cum,tran . These characteristics showed similar variations and trends over the three years. For instance, in 2014, transition day (DOY) from NEE cum source to NEE cum sink (NEE cum,tran ) at ZF was DOY 128, corresponding to CGDD Ta at 104.5 • C. At this point, T s04 was positive but stayed close to 0 • C. This implies that T a plays an important role in NEE cum transition time, as we observed in Figure 2. Before surface T s04 stays nearly 0 • C, the forest was net CO 2 source. Rapid increasing phase in NEE cum followed a steep and a linear increase of CGDD T04 > 4 • C. CGDD T04 increased more than 7 • C in two days (DOY 135-137). After DOY 136, the slope of NEE cum /DOY was −1.05 g C m −2 d −2 . From this point, ZF entered the growing season, and photosynthesis had fully recovered from the cold season. Both photosynthesis and ecosystem respiration continuously increase over the growing season. At ZB, NEE cum,tran was 10 days later (DOY 138) than the NEE cum,tran at ZF. At this time, peats were already warmed up, showing positive temperatures (CGDD Ta = 164.6 • C). Similar to forest, daily mean T s04 also increased approximately 7 • C in two days (DOY 128-130) immediately after CGDD T04 exceeded 1 • C. Presumably, the ground surface at ZF is still partly covered by snow, and shading by trees may result in warmer surface conditions later than at ZB. If mosses were exposed to the open space without shading, then they could receive more direct solar radiation on the surface, which may result in earlier soil warm-up conditions. Characteristics of NEE cum,tran and heat accumulation explained by CGDD Ta and CGDD T04 imply that the bog seems to have a higher level of temperature accumulation than forest, because soil microbial activity requires a suitable temperature level, corresponding to daily T s04 >~7 • C and CGDD T04 of~41 • C. In 2015, NEE cum,tran at ZF occurred in DOY 134, corresponding to CGDD Ta at 137.6 • C. Although the spring mean T a in 2015 was the highest, NEE cum,tran in 2015 occurred the latest compared with other years. Perhaps, warm temperature conditions increase ecosystem respiration, leading to the forest remaining a net CO 2 source for longer. To find a reasonable answer, we will intend to examine components of the flux partitioning for a future study. Anomalously high T a also led to the highest T s04 . CGDD Ta during DOY 100-150 matched with the transition point of CGDD T04 >~4 • C. Compared with ZF, CGDD T04 corresponding to NEE cum,tran at ZB varied 40-75 • C. CGDD T04 at ZB showed more rapid enhancement in other years, possibly associated with the soil microbial responses to warm temperature conditions.

Identifying the Major Drivers of CO 2 Fluxes
At our sites, spring CO 2 fluxes were mainly controlled by air temperature (Table 2). A rapid rise in CO 2 uptake occurred after surface snow melted, and T s04 warmed above 1 • C (Figure 2). We confirmed that abiotic controlling factors of CO 2 fluxes identified from the MARS model agreed with previous findings. In boreal coniferous forests, spring photosynthesis recovery is related to air temperature, surface soil temperatures, first day of snowmelt, and final day of snowmelt [20,[22][23][24][25][26]29]. The coniferous forest was already a net CO 2 sink before the end of snowmelt and while part of the soil surface was still frozen (Figure 2). In a wider context, this result agrees with the finding of Parazoo et al. [55], who showed that spring thaw was the crucial trigger for the start of spring photosynthesis and net CO 2 uptake in boreal forests in Arctic ecosystem. In northern peatland, peat temperature and snowmelt completion have been shown to affect the rapid increase in vegetation productivity [11,18,21,25]. In addition, PAR and temperature have been identified as controlling factors on the CO 2 flux of peat during the growing season [26,56,57].
In our study, T s04 was identified as the primary driver of CO 2 fluxes for both ecosystems (Table 2). From a biophysical point of view, light (i.e., PAR) is known to be the primary driver of the beginning of net CO 2 sink or photosynthesis. In a previous study, a datadriven approach similar to the one used here identified PAR as the most important driver of CO 2 fluxes at the ZF site [42]. During the growing season, including the snow-melting period, PAR strongly controls the growth of Sphagnum [57]. In the boreal winter-to-spring transition period, snow cover decreases as temperature rises, as shown in reverse patterns of surface reflectance and air temperature [21][22][23]. This reflects that, in reality, T s04 and Alb are indirect drivers of CO 2 fluxes.
We chose the MARS model to identify the relative importance of controlling factors of CO 2 fluxes because of its simplicity. In our study site, we often met a long-term data gap, particularly in the winter and early season, meaning that estimates of a reliable annual carbon budget with uncertainty are challenging. To overcome this limitation, further efforts for utilizing major variables identified from the MARS model and/or incorporating with other machine-learning methods suitable for a relatively smaller dataset (e.g., support vector machine) are desirable. With such efforts, we can develop an advanced version for the site-specific gap-filling algorithm.

Potential Drivers of Spring CO 2 Fluxes
We observed that T a , T s04 , Alb, and PAR affect seasonal variations of CO 2 fluxes at a boreal forest and a bog. However, there is still room for exploring unknown or missing drivers of spring CO 2 fluxes. For instance, a recent study by Koebsch et al. [58] suggested that biological drivers seem more important for determining the variability in maximum gross primary productivity compared with abiotic drivers in boreal peatlands. Peichl et al. [59] also emphasized the important role of phenology on seasonal variation of photosynthesis in European peatlands. In addition, it is likely that the rapid increase in net CO 2 uptake after the snowmelt completion at ZB (Figure 2) may be associated with increases in peat temperatures in the rooting zone (0.1-0.2 m), leaf nitrogen, and chlorophyll a concentration [25]. This suggests that chlorophyll a could be a potential biochemical driver of spring CO 2 fluxes in peatlands.
In peatland, hydrological drivers (e.g., precipitation, soil moisture, and water-table depth) and their regime are also crucial for understanding photosynthesis-related processes. Previous studies [11,60,61] suggest that the water-table depth is an important control on the growing season NEE and the annual CO 2 balance, although this control appears to be site dependent. In contrast to previous studies, Strachan et al. [62] did not find significant relationships between peatland productivity or cumulative CO 2 exchange and early season temperature, the timing of snowmelt, or growing season length at an ombrotrophic bog located in Canada. Our data was limited in analyzing the roles of hydrological factors on seasonality of CO 2 uptake in peatland; thus further efforts to examine the roles of biotic and abiotic drivers in regulating peatland C cycle at ZB are desirable. In addition, measurement of other important components of C flux, such as methane, are necessary to better understand the annual C balance in peatland [1,17,63,64].

Role of Air and Soil Temperatures on Vegetation CO 2 Uptake Capacity in Northern Eurasia
Consistent with characteristics of abiotic drivers of CO 2 fluxes shown in Figure 2, NEE cum,tran at ZF requires suitable air temperature accumulation (CGDD Ta >~100 • C) more than soil temperature (Figure 4, Table 5). In previous studies, 5-day running mean of daily mean air temperature and cumulative temperature were good predictors of the commencement of spring photosynthesis in boreal coniferous forests [22,23].
We confirmed that both ecosystems have a zero-curtain period at T s04 that remained close to 0 • C around until mid-April ( Figure 2). This feature is a typical phenomenon in high-northern latitude or alpine ecosystems during the winter-to-spring transition period [21,23,25,60,[65][66][67]. Although the surface snow may not have fully disappeared, warm spells may affect the timing of the start of net CO 2 uptake in boreal forests ( Figure 2). Such phenomena may be a piece of evidence that Scots pines growing in high latitudes have adapted to water-limited environments in severe and long winters. Generally, low soil-temperature conditions during spring constrain water uptake and root activity in boreal conifers [24,66,68,69]. However, there is observational evidence that pine trees in winter use stem-stored water regardless of snowmelt termination or available soil water to maintain their metabolism [22,70]. It is also possible that different tree species may have different strategies for the process of spring photosynthesis recovery [71].
Cold weather appears to temporarily reduce CO 2 uptake rates over the study period. However, plants recovered CO 2 uptake rates as soon as T a increased above 0 • C (Figure 2). We found a similar pattern in other years (Supplementary Figure S1), which tend to be more distinct in the forest than in the bog. In addition, our data reveal that colder spring weather conditions resulted in a reduction of the maximum photosynthetic capacity (Figure 3). In previous studies, intermittent cold temperatures, for instance, frost events, delayed the photosynthesis recovery process of conifers, and the effect is more pronounced if frost is severe [20,72,73]. In peatlands, chilling can reduce chlorophyll a concentration [25,74]. It is unclear to find these features from daily CO 2 fluxes in our study; however, Figure 3 and Table 4 seem to show clues; cold weather temporarily suppresses photosynthetic capacity during spring. Overall, our data showed reasonable ranges and PAR-NEP fit parameters compared with the previous study [64].
A warmer climate may result in an earlier start of spring photosynthesis and increased vegetative productivity as well as annual CO 2 uptake in boreal coniferous forests [29,30,73]. However, the effect of warm temperatures on the net CO 2 uptake capacity of peatlands is likely to be more complex than that of forests because of soil hydrological conditions. Precipitation, water table depth, vegetation type, and phenology are recognised as other important controls on CO 2 fluxes, and their effects on net CO 2 uptake seem to be highly site dependent [11,12,14,37,59].
In this study, both ecosystems showed the highest vegetation productivity in 2015 ( Figure 5). Despite the small amount of data and highly deviated data in 2015, both ecosystems showed negative relationships between NEE cum and mean T a in May. Based on the spatial distribution of 2m temperature anomalies over the study period (2013-2017), it is highly likely that vegetations in Zotino absorb the highest CO 2 uptake among the years because of anomalously warm temperature in May (Supplementary Figure S2). Anomalous warm temperatures were spread over the western Siberia, including the Zotino site. A study by Alekseychik et al. [37] seems to support that our finding is possible. They reported that the spring in 2015 was warmer; monthly mean air temperature in May was 4.1 • C higher than the long-term average in a similar type of bog in Finland. Besides, Liu et al. [14] and Pulliainen et al. [29] found that net CO 2 uptake in boreal ecosystem during the warm spring was greater than the normal year, implying that vegetation is sensitive to the changes in air temperature, and it is a crucial driver to understand spring C uptake. To obtain the robust relationship between NEE cum and mean T a in May, we plan to analyze the longer term of vegetation indices using remote sensing data together with flux data. In addition, Alekseychik et al. [37] also found that despite the very warm spring in 2015, cloudy and rainy weather conditions in summer in that year resulted in a lower CO 2 uptake for the growing season compared with other years. This implies that the high vegetation productivity induced by the anomalously warm temperatures in May 2015 is likely to have been compensated by the cool and wet ensuing summer. Our forthcoming research will examine the annual net CO 2 balance through comparison of the results of the present study with those from other sites in Siberia and for different weather conditions.

Conclusions
CO 2 flux observations at the Zotino sites showed that distinct differences exist in flux magnitude, the timing of start of net CO 2 uptake, and the rate of CO 2 uptake between the boreal forest and bog. In spring, the boreal forest generally changed from net CO 2 source to net CO 2 sink approximately 1-2 weeks earlier than the bog. We found similar aspects for the transition from NEE cum source to NEE cum sink. Rapid net CO 2 uptakes for both ecosystems covaried with T a , corresponding to T s04 values > 5 • C. To change into complete NEEcum sink, ZF seems to require CGDDTa of~80 to 137 • C and ZB requires CGDDTa of 141 to 211 • C. Compared with ZF, seasonal variations in NEE-cum at ZB were smaller. This implies that forests appear more sensitive to the changes in air temperature than bogs under the same spring weather conditions. We confirmed that cold spring weather reduced the maximum photosynthetic capacity in both ecosystems. During the study period, abnormally warm air temperatures in spring 2015 resulted in the highest vegetation productivity. Overall, results suggest that the CO 2 uptake capacities of boreal forest and bog are sensitive to rising air temperature in springtime. This strong relationship between air temperature and NEE seems to support that a linear relationship between NEE estimated from the atmospheric inversion method and air temperature proposed by Rödenbeck et al. [75].
Our analysis is limited in analyzing NEE; thus, future work will examine the effects of anomalously warm spring weather on both photosynthesis and respiration with respect to the annual C balance. Besides, we will investigate the linkage between hydrological conditions (e.g., snow depth and rainfall) in the previous winter and seasonal/annual CO 2 balances in the following year. To overcome the limited number of EC flux measurement, remote-sensing-based phenology and photosynthesis products will be utilized. We anticipate direct CO 2 flux measurements at Zotino will be valuable for evaluating biogeochemical models. Further, it will provide insights into the prediction for the future terrestrial C cycle in northern Eurasia, particularly in remote and relatively undisturbed natural ecosystems.