Thermal acclimation of plant photosynthesis and autotrophic respiration in a northern peatland

Peatlands contain one-third of global soil carbon (C), but the responses of peatland ecosystems to long-term warming are not well understood. Here, we pursue an emergent understanding of warming effects on ecosystem C fluxes at peatlands by constraining a process-oriented model, the terrestrial ECOsystem model, with observational data from a long-term warming experiment at the Spruce and Peatland Responses Under Changing Environments site. Model-based assessments show that ecosystem-level photosynthesis and autotrophic respiration exhibited significant thermal acclimation, with temperature sensitivities being linearly decreased with warming. Using the thermal-acclimated parameter values, simulated gross primary production, net primary production, and plant autotrophic respiration (R a), were all lower than those simulated with non-thermal acclimated parameter values. In contrast, ecosystem respiration simulated with thermal acclimated parameter values was higher than that simulated with non-thermal acclimated parameter values. Net ecosystem CO2 exchange was much higher after constraining model parameters with observational data from the warming treatments, releasing C at a rate of 28.3 g C m−2 yr−1 °C−1. Our data-model integration study suggests that peatlands are likely to release more C than previously estimated. Earth system models may overestimate C uptake by peatlands under warming if physiological thermal acclimation of plants is not incorporated. Thus, it is critical to consider the long-term physiological thermal acclimation of plants in the models to better predict global C dynamics under future climate and their feedback to climate change.


Introduction
Peatlands play an important role in the global carbon cycle.Peatlands are a major global carbon (C) sink, containing one-third of global soil C due to the slow decomposition of organic matter under cold, water-saturated, and oxygen-limited conditions (Bridgham et al 2006).On the other hand, peatlands emit a substantial amount of methane (CH 4 ) into the atmosphere each year (Abdalla et al 2016, Ma et al 2021).The exclusion of peatland C dynamics, including shifts in greenhouse gas emissions, from the majority of Earth system models (ESMs) involved in the Coupled Model Intercomparison Project Phase 6 may contribute to the large uncertainty in ESM projections and low confidence in the magnitude of global soil C losses under warming (IPCC 2021).Therefore, it is critical to understand how peatland ecosystem processes respond to long-term (years to decades) warming so that we can better predict future global carbon fluxes between the atmosphere and Earth's surface and their feedbacks (Smith and Dukes 2013).
Some individual ecological processes of peatlands under long-term warming, e.g.fine-root growth (Malhotra et al 2020), nutrient cycling (Iversen et al 2022), and incorporation and degradation of plant-and microbe-derived organic matter (Ofiti et al 2022), have recently been studied owing to a long-term, whole-ecosystem climate change experiment, i.e. the Spruce and Peatland Responses Under Changing Environments (SPRUCEs) experiment in a peatland ecosystem in Minnesota, USA.Carbon stocks and fluxes at the ecosystem level, including methane emissions, in response to long-term warming have also been investigated under this project (Griffiths et al 2017, Ma et al 2017, Jiang et al 2018, Hanson et al 2020, Ricciuto et al 2021, Yuan et al 2021).Net ecosystem production shifted significantly toward a smaller C sink or a minor C source under the most extreme experimental warming scenario (i.e.+9 • C warming scenario; Hanson et al 2020).However, peatland ecosystem-level plant physiological responses to long-term warming remain unclear.
Physiological mechanisms for instantaneous temperature responses of C exchange processes of plants are relatively well-studied at the leaf level.For example, it has been widely reported that increased temperature stimulates enzyme activities up to an optimum temperature, which then declines due to enzymatic degradation at higher temperatures (Niu et al 2012, Smith andDukes 2013).Compared to the leaf level, less is known about ecosystem-level C exchange responses to warming, especially in the long term (Bradford et al 2008, Ziehn et al 2011, Lombardozzi et al 2015, Smith et al 2015).Empirical evidence indicates that the initial increase in ecosystem fluxes such as soil or ecosystem respiration (ER) under warming can decline or disappear after a certain time (Luo et al 2001, Melillo et al 2002, Eliasson et al 2005).Such changes in C cycling responses to warming can be attributed to physiological acclimation, species composition change, or altered substrate availability (Bradford et al 2008, Smith andDukes 2013).
The down-regulation of temperature-adapted physiological behavior, which is defined as thermal acclimation, may lead to the observed C flux responses to temperature (apparent acclimation) or slow down the increasing trend of C fluxes (non-apparent acclimation) in a warmer climate (Oechel et al 2000, Davidson andJanssens 2006).Thermal acclimation of plant photosynthesis and autotrophic respiration has frequently been reported (e.g.Campbell et al 2007, Smith and Dukes 2013, Smith et al 2015, 2019, Kumarathunge et al 2019).As a result of such acclimation, simulated C sensitivity to climate was significantly reduced (Smith et al 2015, Liang et al 2018).Therefore, failure to account for acclimation of plant photosynthesis and autotrophic respiration may cause bias in simulating land C cycling (Smith and Dukes 2013, Smith et al 2015, Liang et al 2018, 2019).More and more models have incorporated dynamic parameterization that allows models to simulate the thermal acclimation of plant photosynthesis and autotrophic respiration (Atkin et al 2014, Lombardozzi et al 2015, Huntingford et al 2017).
Photosynthetic thermal acclimation is primarily controlled by photosynthetic biochemistry, although stomatal regulation and daytime respiration also contribute (Lin et al 2013, Kumarathunge et al 2019).The algorithms used for approximating photosynthetic acclimation to temperature emerged only recently, primarily due to the lack of empirical data for temperature-responsive variables such as the maximum carboxylation rate at reference temperature (V cmaxref ), usually 25 • C and the maximum potential rate of electron transport at a reference temperature (J maxref ), which are most commonly used in implementations of the Farquhar photosynthesis model (Farquhar et al 1980) in most ESMs (Kattge and Knorr 2007, Friend 2010, Scafaro et al 2017).V cmaxref and J maxref represent basal values of photosynthetic capacity while activation energy (E a ), deactivation energy (E b ), and entropy (∆S) describe the direct effect of temperature on photosynthetic capacity (Stinziano et al 2018).Hikosaka (2006) identified two responsive variables for photosynthetic thermal acclimation of 23 C 3 species based on empirical data: activation energy of V cmaxref (E aV ) and J maxref (E aJ ), as well as the ratio of J maxref to V cmaxref (J/V).Later, linear regressions of J/V and entropy of V cmax (∆S V ) and J max (∆S J ) against plant growth temperature (T growth ) were derived for 36 species, whereas no significant correlation between E a and T growth among these species was found (Kattge and Knorr 2007).The generalized linear regression models have widely been used in large-scale and global modeling studies (Ziehn et al 2011, Arneth et al 2012, Lombardozzi et al 2015, Smith et al 2015).The most comprehensive analysis of the photosynthetic temperature responses covered 141 species from tundra to tropical forest (Kumarathunge et al 2019), in which ∆S V , ∆S J , J maxref , and J/V all declined but V cmaxref did not change with increasing T growth .Kumarathunge et al (2019) also found an increase in E aV but no change in E aJ with increasing T growth .The aforementioned studies described thermal acclimation as a function of T growth and responsive variables, such as J/V and temperature sensitivity (Q 10 ).More empirical data are needed to generalize the key model parameters (e.g.V cmaxref ) to accurately simulate long-term feedbacks between the atmosphere and terrestrial ecosystems (TECOs) under climate warming.
Dynamic parameterization has also been incorporated to allow models to simulate the thermal acclimation of plant autotrophic respiration (Atkin et al 2014, Lombardozzi et al 2015, Huntingford et al 2017).Early studies aimed to represent the dynamic responses of autotrophic respiration to temperature using a temperature-dependent Q 10 , instead of a fixed value, to capture the instantaneous responses of autotrophic respiration observed in the field (Larcher 1980).Other studies used algorithms to characterize respiratory acclimation to longer-term growth conditions (King et al 2006, Atkin et al 2008).The algorithms used by King et al (2006) and Atkin et al (2008) were both derived from empirical data in warming experiments with multiple species that covered different functional groups.However, more empirical research is needed to quantify the dynamic responses of plant respiration to temperature across a wider range of species, plant functional types, and environmental conditions.
The lack of empirical data for identifying the key variables of thermal acclimation has become a barrier to the realistic representations of thermal acclimation of plant photosynthesis and autotrophic respiration in models.While more observations and experiments are needed for deriving the temperature-responsive variables, an alternative approach, data-model fusion, can help obtain those temperature-responsive variables.Data-model fusion assimilates empirical data from multiple sources into process-based biogeochemical models to constrain model parameters (Luo and Schuur 2020), including parameters describing thermal acclimation of plant physiological processes (e.g.basal photosynthetic rate and Q 10 ).For example, a data-model fusion study in Arctic tundra showed that warming increased the light-use efficiency of vegetation but decreased the baseline turnover rate of both labile and recalcitrant soil organic carbon (SOC) pools, suggesting that physiological acclimation in plants and functional gene shifts in microbes occurred (Liang et al 2018).
To examine how ecosystem fluxes of peatlands responded to long-term warming, in this study, we assimilated observational data from the SPRUCE experiment to a process-based biogeochemistry model, the TECO model, specialized for the SPRUCE site (TECO_SPRUCE).We hypothesized that plant physiological processes such as photosynthesis and autotrophic respiration exhibit thermal acclimation under long-term warming.

Site description and measurements
The empirical data used in this modeling study were measured from the SPRUCE experiment (Hanson et al 2017; http://mnspruce.ornl.gov).The experiment site is a precipitation-fed, ombrotrophic bog in northern Minnesota, USA (N 47 • 30.476 ′ , W 93 • 27.162 ′ ).The mean annual temperature was 3.4 • C and the mean annual precipitation was 780 mm from 1961 to 2009 (Sebestyen et al 2011).The mean peat depth is 2-3 m due to historical accumulation of peat under cold, anaerobic conditions (Parsekian et al 2012).The site includes the dominant tree species, Picea mariana and Larix laricina, a variety of ericaceous shrubs, such as Rhododendron groenlandicum and Chamaedaphne calyculata, herbaceous perennials Eriophorum vaginatum, Carex trisperma, and Maianthemum trifolium, as well as a dense layer of Sphagnum sp.moss.The herbaceous plants have seasonal dieback of their aboveground tissues.
An ecosystem-level manipulation was implemented at this site to study the responses of the northern peatland ecosystems to climate warming and elevated atmospheric CO 2 (Hanson et al 2017).For this specific study, we focused on warming treatments only.
A gradient of warming treatments, including controls, ambient (+0), +2.25, +4.5, +6.75, and +9 • C warming is implemented in open-top chambers (12 m in diameter and ∼8 m in height).The control refers to the plot with no open-top chamber whereas the ambient refers to the plot covered by an open-top chamber but with no warming.Each of the warming enclosures is heated both above-and below-ground at the same target temperatures.Deep peat warming began in June of 2014, and target temperatures were achieved at a depth of 2 m.The aboveground warming started in August of 2015, and the target temperatures were controlled at 2 m above ground.Air temperature in the control plots without chambers is approximately 2 • C lower than the ambient plots with chambers but no heating due to the heating effect of the chamber.Thus, the ambient plots are considered the lowest warming treatment in this study and are named +0 • C. The environmental variables measured in 2011-2014 from the earliest control plot were used to spin up the TECO_SPRUCE model to a steady state.The environmental variables from each warming plot from 2015 to 2018 were used as model input (table 1).Instead of simulating temperature in different peat layers, we used the measured temperatures half-hourly at 0, −5, −10, −20, −30, −40, −50, −100, and −200 cm depths and aggregated them into hourly intervals to drive the model because the peat was heated to the target temperatures at a 2 m depth.Soil moisture of different soil layers and water table depths from individual plots in 2014-2018 were used to calibrate the modeled water-heat balance.We used ecosystem-level C pools (leaf, wood, root, and SOC) and fluxes, including gross primary production (GPP), net ecosystem CO 2 exchange (NEE), and ER, of the shrub ecosystems from 2015 to 2018 for data-model fusion (Hanson et al 2016a, Norby et al 2019).These fluxes were measured at least once a month during snow-free months from 2015 to 2018 with a portable open-path analyzer attached to a chamber (1.2 m in diameter).The chamber was put on the soil surface within the plots and covered the shrubs.Data points measured within an hour were averaged as hourly measurements for each plot to match model's time step for data-model fusion as well as model evaluation.Finally, about 30 data points from each plot were used to compare with hourly model output for data-model fusion.A complete list of the data streams used as observations to constrain model parameters is in table 1.

Model description
We used a process-based biogeochemistry model, TECO_SPRUCE, for our data-model fusion study.The model was built with six major modules working at an hourly time step: canopy photosynthesis, soil water dynamics, plant growth, soil thermal dynamics, soil carbon/nitrogen (N) transfer, and soil CH 4 dynamics.A detailed description of these modules can be found in Weng and Luo (2008), Shi et al (2015), Jiang et al (2018), Huang et al (2017) and Ma et al (2017Ma et al ( , 2022)).Here, we briefly describe these modules and highlight the instantaneous temperature response functions for photosynthesis and autotrophic respiration.
The canopy photosynthesis module was mainly derived from a two-leaf model which coupled surface energy, water, and carbon fluxes (Wang and Leuning 1998).Leaf photosynthesis was calculated based on the Farquhar photosynthesis model (Farquhar et al 1980) and the stomatal conductance model (Ball et al 1987).The temperature responses of parameters V cmax (maximum rate of carboxylation) and J max (maximum potential rate of electron transport) varied with leaf temperature following the modified Arrhenius function (Johnson et al 1942): (1) where V cmaxref is V cmax at reference temperature (293.2K), T leaf is leaf temperature (K), f (T leaf ) is the factor of leaf temperature (unitless), J/V is the ratio of J maxref to V cmaxref , ∆H a is the activation energy (J mol −1 ), T ref is reference temperature of leaf (K), R gas is the gas constant (8.314J K −1 mol −1 ), ∆S is entropy (J K −1 mol −1 ), and ∆H d is the deactivation energy (J mol −1 ).∆S used in equations ( 1) and ( 2) have slightly different values, where Entropy of J max (∆S j ) is derived from Entropy of V cmax (∆S v ): ∆S j = ∆S v * 668/664.According to Kattge and Knorr (2007), the photosynthetic acclimation to growth temperature is reflected by changing values of ∆S and J/V.Thus, in this study, instead of using a constant value for ∆S and J/V, we set these parameters to be variable within a prior range, which could be constrained by the observation data.The water table level was estimated using a simple bucket model as described by Granberg et al (1999).The plant growth module calculated the allocation of photosynthesis C to plant pools (leaf, stem, and root), plant growth, plant autotrophic respiration, phenology, and C transfer to litter and soil pools.The autotrophic respiration of leaf, stem, and root was calculated as: where Rm leaf , Rm stem , and Rm root are maintenance respiration rates of leaf, stem and root, respectively (g C m −2 h −1 ), R 0leaf , R 0stem , and R 0root are basal respiration rates of leaf, stem and root, respectively (g respired C g −1 biomass C m −2 h −1 ), S NRa is the nitrogen scaler for autotrophic respiration (unitless), leafC, stemsapC, rootsapC are C content of leaf, stem sapwood, and root sapwood, respectively (g C m −2 ), Q 10Ra is temperature sensitivity of autotrophic respiration (unitless), Tair is air temperature, and fnsc is the scaling factor of nonstructural C pool (unitless).We set prior ranges for basal respiration rates of leaf, stem, and root and temperature sensitivity of autotrophic respiration and they could be constrained by the observation data.

Data-model fusion
We applied the Metropolis-Hasting algorithm (Metropolis et al 1953) to generate the posterior distribution of parameters.We assumed a uniform distribution for the prior range of each parameter so that the chance of each value being accepted is equal.We also assumed that errors between observations and model simulations independently follow a normal distribution with a zero mean.The cost function weights the mismatch between the multiple observational data streams and the modeled corresponding variables, represented by: where Z i (t) is the ith observation stream at time t, X(t) is the corresponding simulated value, and σ i (t) is the standard deviation of observational error estimates.There were seven observational data streams used in total (table 1).We generated the posterior distributions of parameter values using 50 000 iterations during the optimization process.The parameter value at the current step was based on the accepted parameter value in the previous step.The current value was accepted only if the mismatch between the observation and model simulation was reduced or otherwise randomly accepted with a 0.05 probability.We used the Gelman-Rubin statistic (Gelman and Rubin 1992) to examine the convergence of sampling chains.The first half of the accepted parameter values were discarded from the burn-in period, and the second half of the accepted parameter values were used for posterior analysis.More details on sampling and the cost function can be found in Xu et al (2006).
Under photosynthetic thermal acclimation, the ratio of J maxref and V cmaxref (i.e.J/V), ∆S V , and ∆S J have consistently been found to change with growth temperature across different studies, but activation energy (∆H a ) did not change with growth temperature (Hikosaka et al 2006, Kattge and Knorr 2007, Kumarathunge et al 2019).Given the limited amount of observational data to constrain the model parameters, if we allowed all the six parameters (J maxref , V cmaxref , ∆S V , ∆S J , ∆H aV , and ∆H aJ ) in equation ( 7) to vary during optimization, the parameter posterior distributions tend to be highly correlated to each other, resulting in poor results.We thus set ∆H aV and ∆H aJ to be constant and let J maxref , V cmaxref , ∆S V , and ∆S J vary.In addition to these four parameters, we selected all other instantaneous temperature-responsive parameters (table 2) that might change due to thermal acclimation according to previous studies to be constrained by observation data.Additionally, parameters related to changes in C pool size, such as turnover time, allocation factor from/to labile soil C pool, and maximum growth rate of leaf, stem, and root C, were also selected for data-model fusion.We selected 15 parameters in total and detailed information on these parameters and their prior ranges are listed in table 2. The prior ranges of parameters were based on empirical measurements or modeling studies from peatland ecosystems.Indeed, in our test runs of the data-model fusion, we also included all the allocation parameters of plant C pools to vary but having too many parameters to vary at the same time might obscure the limited information contained in observation data streams, especially with limited data points.We thus set constant allocation parameters of plant C pools.We randomly saved 500 sets of the simulation results from each data-model fusion and calculated the mean value and standard deviation of each parameter to analyze the responses of ecosystem C fluxes to long-term warming and the potential shifts in the parameter values due to thermal acclimation.We further explored how changes in parameters due to physiological thermal acclimation could influence the ecosystem C fluxes of peatlands.Here we used 'thermal acclimated parameters' and 'non-thermal acclimated parameters' to refer to parameter values constrained by assimilating observation data from warming (indicated by color lines in figure 1) and control plots, respectively.Specifically, we ran the model with both non-thermal acclimated and thermal acclimated parameter values and compared the results of the two runs.Both runs used the initial conditions in the control plot, driven by six sets of drivers corresponding to six warming treatments.

Statistical analyses
To evaluate model performance in simulating C fluxes, linear regressions were performed to compare simulated and observed C fluxes.The observed C fluxes used for evaluation are the same C fluxes as used for constraining the TECO model, which were measured 1-2 times a month during growing seasons from 2015 to 2018 in each plot.Root mean square error (RMSE) was calculated to evaluate model performance against the observation data.
To test the possible thermal acclimation of plant physiology, we performed linear regressions between model parameters and temperature, using yearly average of parameters and temperature, including all control and warming treatments from 2016 to 2018.Lastly, we analyzed the dependency of the simulated annual C fluxes on the mean annual temperature including all treatments from 2016 to 2018, in which actual mean annual air temperatures were used as independent variables.All statistical analyses were performed with R software (version 3.6.1,R Core Team 2017).

Parameters constrained by data-model fusion
Among the 15 parameters, eight parameters, J/V, V cmaxref , ∆S V , ∆S J , Q 10Ra , gddonset, τ leaf , and τ root , were well constrained by the observational data streams across all treatments, with a unimodal-shaped posterior distribution (table 2, figure 1).The model using optimized parameters by data assimilation fitted the observational data well with RMSE between 0.047 and 0.126 (figure 2).
The correlations between parameter mean values of the posterior distributions and mean annual temperature including all treatments from 2016 to 2018, during which the whole-ecosystem warming treatments were active all year round, were shown in figure 3. The instantaneous temperature responsive parameters of photosynthesis, including V cmaxref , J maxref (which was calculated by multiplying J/V by V cmaxref ), entropy of V cmaxref (∆S V ), and entropy of J maxref (∆S J ) all linearly decreased with the increasing treatment temperature.Despite being well constrained, J/V did not show a significant correlation with the increasing mean annual air temperature.
While the basal autotrophic respiration rates of leaf, stem, and root were not well constrained (table 2), the temperature sensitivity of autotrophic respiration, Q 10Ra , was well constrained and the mean values of the posterior distribution across all treatments exhibited a significant linear decrease with increasing mean annual air temperatures at a rate of −0.047 • C −1 (P < 0.001, figure 3).

Simulated ecosystem carbon fluxes under warming
The comparison of C fluxes from 2016 to 2018 simulated with non-thermal acclimated parameter values vs. thermal acclimated parameter values were shown in figure 4. Simulated GPP (figure 4(a)) and net primary production (NPP) (figure 4(b)) using non-thermal acclimated parameter values were higher than those using thermal acclimated parameter values.Moreover, GPP simulated with thermal acclimated parameter values did not change significantly with increasing annual air temperature.In contrast, NPP simulated with thermal acclimated parameter values increased linearly with increasing annual air temperature at a rate of  19.0 g C m −2 yr −1 • C −1 (P < 0.001).Similar to GPP and NPP, simulated plant autotrophic respiration, R a , using non-thermal acclimated parameter values was higher than that using thermal acclimated parameter values (figure 4(c)).R a simulated using thermal acclimated parameter values decreased linearly as annual air temperature increased by 17.2 g C m −2 yr −1 • C −1 (P < 0.001).
Simulated ER with both non-thermal acclimated and thermal acclimated parameter values increased significantly with annual air temperature, and ER simulated with thermal acclimated parameter values was higher than ER simulated with non-thermal acclimated parameter values (figure 4(d)).Without thermal acclimation (simulations with non-thermal acclimated parameter values), due to significant increases in both GPP and ER at a similar rate, the ecosystem remained a minor C sink, i.e. a negative value of NEE (figure 4(e)).However, after constraining parameters with observational data from the warming treatments, the ecosystem turned from a neutral/minor C sink into a major C source, releasing C at a rate of 28.3 g C m −2 yr −1 • C −1 (P < 0.001; figure 4(e)).

Discussion
Peatlands are an important ecosystem in the global C cycle as they store large amounts of soil C, which could be vulnerable to climate change.However, the responses of peatland C cycle to climate change, especially long-term responsive patterns, are poorly understood.In this study, leveraging measured data from a unique manipulative peatland experiment with multiple warming treatments to constrain a process-based ecosystem model, we explored ecosystem C flux response to long-term warming.We found that including representing thermal acclimation in model processes significantly altered ecosystem C flux components in response to warming; ultimately leading to the ecosystem becoming a much higher net C source than if thermal acclimation was excluded from the model.Below we discuss our findings and their implications for land C modeling.

Thermal acclimation of plant photosynthesis and autotrophic respiration at ecosystem level
Significant linear responses of plant photosynthesis and autotrophic respiration variables, such as J/V and entropy, to increasing growth temperature have been found at the species level (Kattge and Knorr 2007).These linear temperature responses of photosynthesis are affected by a variety of factors (e.g.soil moisture, slow vs. fast-growing strategies, and plant functional types), which are not always been taken into account when reviewing photosynthetic acclimation to temperature across studies (Lin et al 2013).Similarly, the ability of plants to adjust their respiratory rates in response to long-term thermal changes also varied by location (Larigauderie and Körner 1995) and plant functional types (Tjoelker et al 1999, Atkin et al 2006).The linear responsive pattern drawn from the existing studies supports the theory that plants maximize their resources for growth and reproduction through acclimation, and climate alone can be used to predict the optimal rate of photosynthesis and autotrophic respiration of plants (Smith et al 2019).However, such linearity has never been directly tested at the ecosystem level at one site with a gradient of warming treatments.The unique design of the SPRUCE experiment and the use of the data-model fusion approach makes it possible to test whether model parameters related to C fluxes responded linearly or nonlinearly to increasing temperature at the ecosystem scale.We found that after four years of warming treatments, photosynthetic parameters sensitive to temperature, including V cmaxref , J maxref , and entropy, exhibited thermal acclimation, decreasing linearly with increased temperatures (figure 3).In the meta-analysis by Kattge and Knorr (2007) based on 36 different C 3 species from various ecosystem types, they found that the entropy of both V cmaxref (∆S V ) and J maxref (∆S J ) declined as growth temperature increased, at a rate of 1.07 and 0.77 J K −1 mol −1 • C −1 , respectively.A more comprehensive study based on a wide distribution of 141 species, from tundra to tropical forest, concluded that both ∆S V and ∆S J decreased linearly with an increase in temperature by 1.5 J K −1 mol −1 • C −1 (Kumarathunge et al 2019).Our results based on the regressions between mean values of the posterior distributions against mean annual temperature also showed a significant linear reduction in ∆S V with a slope of −1.18 J K −1 mol −1 • C −1 and in ∆S J with a slope of −1.19 J K −1 mol −1 • C −1 , both of which lie in the middle of the empirical-based values.
Unlike previous studies, we did not find a significant change in the ratio of J maxref and V cmaxref in response to the change in temperature.Instead, we did find a linear decrease with increased mean annual temperature in both J maxref (−0.76 µmol m −2 s −1 • C −1 ) and V cmaxref (−0.48 µmol m −2 s −1 • C −1 ).The different results for the simulated responses of photosynthesis to temperature change in our study compared with previous studies could be attributed to the adoption of a decreased J/V at the higher temperature, because this ratio is a mathematical term, but it is J maxref and V cmaxref that directly determine the photosynthetic rate.Reducing J maxref alone is likely to reduce total photosynthesis (Arneth et al 2012, Lombardozzi et al 2015, Smith et al 2015), whereas increasing V cmaxref alone is likely to increase total photosynthesis (Lin et al 2012).According to Mercado et al (2018), a simultaneous decrease in both J maxref and V cmaxref may be the best way to achieve a decreased J/V.Although we did not find a significant decrease in J/V, our results showed that both J maxref and V cmaxref declined linearly as mean annual temperature increased, which has been proposed by Mercado et al (2018) to account for photosynthetic thermal acclimation.
In our study, the temperature sensitivity of plant autotrophic respiration, Q 10Ra , was well-constrained across all the treatments and the mean values of the posterior distribution of Q 10Ra declined linearly with an increase of mean annual air temperature at a rate of 0.047 • C −1 (figure 3).Based on results from 56 species that covered arctic, boreal, temperate, and tropical biomes, Tjoelker et al (2001) found a negative linear regression between temperature sensitivity of autotrophic respiration and increased growth temperature, with a slope of −0.046 • C −1 and an intercept of 3.22.Later, another study with 121 species also found a similar relationship with a slope of −0.043 • C −1 and an intercept of 3.09 (Atkin and Tjoelker 2003).Our modeling analysis at the ecosystem level from the northern peatland ecosystem is in line with the findings observed in these studies, and also supports the theory that climate alone can be used to predict the optimal reaction rate, irrespective of plant functional types (Smith et al 2019).Physiological acclimation we found in this study may be associated with acclimation of newly emerged and overwintered leaves, but not mature leaves of shrub species in this ecosystem (Ward et al 2019).

Ecosystem carbon fluxes under long-term warming
Photosynthetic thermal acclimation can strongly affect ecosystem-atmosphere feedbacks in the global C cycle, especially as the climate warms.The incorporation of thermal acclimation of photosynthetic parameters into a canopy flux model can improve modeled GPP when compared with eddy covariance data (Stinziano et al 2018).By fusing the observational C fluxes into the TECO_SPRUCE model, we found the physiological thermal acclimation of V cmaxref , J maxref , and entropy in photosynthesis, which resulted in the apparent thermal acclimation of GPP (figure 4(a)).In comparison, if physiological thermal acclimation was not considered, the model would simulate a much higher GPP in response to temperature increase with an increased rate of 22.5 g C m −2 yr −1 • C −1 (figure 4(a)).In addition to warming, there are a few other covarying factors that might have contributed to the down-regulation of GPP.First, the lower relative humidity in the warmer treatments caused moss desiccation and consequent smaller coverage of moss (Desai 2014, Norby et al 2019).Second, nutrient limitation, though not yet reported at this experiment site, might prohibit a positive effect of warming on GPP (Dusenge et al 2019).However, such down-regulation of GPP may not be detected in short-term warming experiments.For example, Johnson et al (2013) found two years' warming treatment with infrared lamps significantly increased GPP of the drier, hummock plots whereas had no significant effects on GPP of the wetter, lawn plots in a northern Michigan peatland.
Similar to photosynthetic thermal acclimation, apparent autotrophic respiration was significantly lower at higher temperatures due to a reduced Q 10Ra .Opposite to GPP, a slight increase in apparent autotrophic respiration rates would be simulated if thermal acclimation of autotrophic respiration was not taken into account (figure 4(c)).The net effect of an unchanged GPP and decreased R a under warming was a slow increase in NPP, which was lower than that without incorporating physiological thermal acclimation (figure 4(b)).Increased aboveground NPP (ANPP) of shrub and community-level belowground net primary production have been observed along a gradient of increasing warming in a manipulative mesocosm experiment with intact soil monoliths taken from a bog in northern Minnesota (Weltzin et al 2000).However, ANPP of graminoid was found to decrease with increasing warming and ANPP of bryophyte was not affected by warming in the same manipulative experiment (Weltzin et al 2000(Weltzin et al , 2001)).Moreover, the responses of NPP to warming in this modeling analysis are inconsistent with that estimated with different components of NPP measured at the same experiment and in the same time period (i.e. 2016-2018), in which warming-induced decreases in aboveground production by trees and the Sphagnum moss community were offset by warming-induced increases in aboveground production by the shrubs and in the belowground production of fine roots of the woody vascular species (Norby et al 2019, Hanson et al 2020, Malhotra et al 2020, McPartland et al 2020).Moreover, the increase of NPP under warming simulated in this study is contrasted with a modeling analysis by Jensen et al (2019).In their study, by using seasonal, species-specific photosynthesis and respiration parameters measured at this SPRUCE site before the treatments were implemented to parameterize a land surface model, ELM-SPRUCE, ELM-SPRUCE predicted decreased NPP under warming treatments for all four species examined, including two understory shrubs (R. groenlandicum and C. calyculata) and two overstory trees (P.mariana and L. laricina).The different responses of NPP to warming among these studies need further investigation, e.g.different measurement methods and/or model intercomparisons.
In contrast to the down-regulation of plant physiological processes, ER, ER, was much higher after adjusting parameters with observational data than that without adjustments of parameters (figure 4(d)).The relative importance of GPP and ER in determining NEE varied.For example, a meta-analysis found that GPP had a greater contribution to NEE than ER under warming (Niu et al 2012).In another study, Valentine et al (2000) reported that ER dominated the variability in the annual C balance of 15 European forests.In this study in a peatland ecosystem, as a result of the combination of physiological down-regulation under warming and altered parameters for ER by warming, the TECO_SPRUCE model projected a C source of this peatland ecosystem in response to warming at a rate of 28.3 g C m −2 yr −1 • C −1 , in comparison to a C sink at a rate of 3.5 g C m −2 yr −1 • C −1 if the model parameters modified by warming treatments were not incorporated.Estimated ecosystem net C exchange with measured components of net C exchange at the same experiment during the same time period also showed a rapid net carbon loss by 31.3 g C m −2 yr −1 • C −1 (Hanson et al 2020) and the carbon loss rate has been found to increase to 34-35 g C m −2 yr −1 • C −1 after the measurements in 2019-2021 are incorporated (personal communications with Dr Paul Hanson).However, the carbon loss in the measured net C exchange was dominated by warming-enhanced decomposition losses of CO 2 and enhanced net CH 4 production under warming treatments because NPP did not change under warming (Hanson et al 2020), which is consistent with the findings by Johnson et al (2013) that the insignificant warming effects on NEE was primarily dominated by the warming effect on ER rather than GPP in the hummock plots in a northern Michigan peatland.Further studies are needed to explore the slightly different mechanisms for the rapid carbon loss between the simulated and the estimated net C exchange with measurements of its components, i.e. the difference in the contribution of NPP to net C exchange.It should be noted that the initial acceleration in CO 2 -C losses via ER in drying and warming in continental bogs may be replaced by the peatland's original CO 2 -C sink function because the persistent drought and warming can shift vegetation composition, resulting in increased NPP over time as revealed by Munir et al (2015) who compared C fluxes among three sites with contrasting water table levels in a dry continental treed bog.

Implications for incorporating thermal acclimation into ESMs
It has been well acknowledged that model parameters need to be calibrated with changing temperatures to represent biological thermal acclimation (e.g.Kumarathunge et al 2019, Lawrence et al 2019, Luo and Schuur 2020).An increasing number of models have begun to incorporate algorithms with dynamic responsive variables that allow models to simulate thermal acclimation of photosynthesis and autotrophic respiration.It is still challenging to parametrize models because of the lack of empirical data available for parameter optimization to implement thermal acclimation of plant photosynthesis and autotrophic respiration.Our study demonstrated an alternative approach to incorporate physiological thermal acclimation by using a process-based ecosystem C cycle model and in situ data.This novel method facilitates explorations of how ecosystem fluxes respond to long-term warming.Using warming-induced parameter adjustments had a significant effect on simulated ecosystem fluxes, shifting the peatland ecosystem from C sink to C source.As such, ESMs may overestimate C uptake under global warming if thermal acclimation is not implemented in the models (Liang et al 2018).The linear responses of model parameters to changes in temperature in this study are comparable with previous studies and can be implemented into ESMs.Certainly, more ecosystem types need to be investigated in the future to test these algorithms and to collect data needed to better simulate how the global C cycle responds to climate warming and its feedback to climate.

Conclusions
In this study, by assimilating observational data from an existing experiment at SPRUCE site involving a gradient of warming treatments into a process-oriented TECO_SPRUCE model, we examined how ecosystem fluxes in the peatland responded to long-term warming.We found thermal acclimation of photosynthesis and autotrophic respiration of plants at the ecosystem level, demonstrated by a linear decrease in temperature sensitivities with warming.The down-regulation of photosynthesis, in combination with increased ER under warming, led to a shift from a C sink to a C source in this peatland ecosystem.Our results at the ecosystem level were consistent with the findings previously observed at the species level and supported the theory that climate alone is a good predictor of the optimal rate of photosynthesis and autotrophic respiration of plants.Without incorporating the thermal acclimation of photosynthesis and autotrophic respiration of plants, ESMs may overestimate the C sequestration capacity of peatlands under warming.Thus, it is critical to incorporate the long-term thermal acclimation of plant physiological processes into ESMs to better predict global C dynamics under future climate and their feedbacks to climate change.

SFigure 1 .
Figure 1.Posterior distributions of parameters from Metropolis-Hasting simulations by assimilating observation data from control and warming plots.Parameter definitions are listed in table 2. ∆Sv and ∆S j have similar distributions because ∆S j is calculated from ∆Sv multiplied by a constant, see details in equation (3).

Figure 2 .
Figure 2. Comparison of observed and model-simulated (mean ± standard deviation) gross primary production (GPP), net ecosystem CO2 exchange (NEE), and ecosystem respiration (ER) in the control and warming treatments.

Figure 3 .
Figure 3. Parameter mean values from well-constrained posterior distributions to characterize the temperature dependence on mean annual air temperature including all control and warming treatments from 2016 to 2018: (a)-(e) are temperature response variables in photosynthesis, (f) is temperature response variable in autotrophic respiration, and (g)-(h) are the turnover time of root and leaf pools in unit of year.Linear regressions are shown in solid lines.The slope (m) of linear regression indicates the magnitude of the parameter in response to increasing mean annual air temperature at 2 m height.* , * * , and * * * indicate significant levels at 0.05, 0.01, and 0.001, respectively.

Figure 4 .
Figure 4. Simulated annual carbon fluxes in relation to mean annual air temperature including all control and warming treatments from 2016 to 2018.Points are mean values from 500 model simulations with parameters randomly drawn from the posterior distributions.Green and gray points are simulated C fluxes with thermal acclimated and non-thermal acclimated parameter values, respectively.Mean annual air temperatures at 2 m height were used for regressions.

Table 1 .
The observational data measured in the SPRUCE experiment and used in this study.

Table 2 .
Parameters included in data-model fusion.