Temperature-Dependent Oxygen Isotope Fractionation in Plant Cellulose Biosynthesis Revealed by a Global Dataset of Peat Mosses

The oxygen isotope composition (δ18O) of plant cellulose has been widely used to study ecohydrological processes of ecosystems as well as to reconstruct past climate conditions in terrestrial climate archives. These applications are grounded on a key assumption that the biochemical fractionation during cellulose synthesis is a constant around +27‰ and is not affected by environmental factors. Here we revisit the influence of temperature on biochemical fractionation factor during cellulose synthesis using a global compilation of Sphagnum cellulose δ18O data. Sphagnum (peat mosses) are known for inhabiting waterlogged peatlands and possessing unique physiological strategy in that their cellulose δ18O could closely reflect growing-season precipitation δ18O. Although within-site cellulose δ18O variability shows a median standard deviation of 0.7–0.8‰ resulting from different degrees of evaporative enrichment of 18O in metabolic leaf water, this evaporative enrichment should be a small quantity due to the external capillary “water buffer” in Sphagnum mosses. Using site-specific minimum cellulose δ18O data that most likely reflect the signal of unevaporated source water, we show that the apparent enrichment factor between cellulose and precipitation δ18O increases with decreasing air temperature. In particular, the apparent enrichment factor could reach values as high as +32‰ when growth temperature is below 5°C. This observational dataset extends the support for the temperature-dependent oxygen isotope fractionation in plant cellulose synthesis previously demonstrated in laboratory experiment, with implications for paleoclimate and plant physiology studies that employ cellulose δ18O measurements, particularly in alpine regions.


INTRODUCTION
The oxygen isotope ratio ( 18 O/ 16 O, expressed as δ 18 O) of plant cellulose (hereafter δ 18 O c ) is a valuable tracer to characterize the environmental conditions during tissue growth (Barbour, 2007;Sternberg, 2009). It has been used as a proxy to reconstruct past climate conditions in tree rings (Libby et al., 1976), lake sediments (Edwards and McAndrews, 1989), and peat deposits (Hong et al., 2000). Earlier studies established empirical relationships between δ 18 O c and climatic parameters, such as temperature, relative humidity, and precipitation amount for a given region (Libby et al., 1976;Treydte et al., 2006;Managave et al., 2010). An improved understanding of physiological and biochemical processes enabled developing mechanistic models to quantitatively decipher the δ 18 O c signal for paleoclimate reconstruction Kahmen et al., 2011). It has been well recognized that the source water for vascular plants, in most cases precipitation, is taken up by roots for photosynthesis in leaves, while the leaf water is enriched in 18 O relative to precipitation due to plant transpiration through stomata. The carbohydrate product (e.g., sucrose) that has been further enriched in 18 O by biochemical fractionation is then translocated to the site of cellulose synthesis, such as stem, where additional biochemical exchange occurs with stem water. As such, plant δ 18 O c signal essentially reflects precipitation δ 18 O (hereafter δ 18 O p ) with an offset governed by both physiological and biochemical processes Sternberg, 2009).
The current mechanistic models have made progress in depicting the physiological effects including how different climatic parameters affect the degree of 18 O enrichment in leaf water and imprint the δ 18 O c signal (Kahmen et al., 2011). However, less is known about the biochemical effects that are equally important in shaping the final δ 18 O c signal. The combined physiological and biochemical effects could be described by the following equation (Barbour and Farquhar, 2000): where δ 18 O c and δ 18 O p are the oxygen isotope composition of cellulose and precipitation, respectively; leaf is the degree of 18 O enrichment in leaf water known as the physiological effect; p x is the proportion of unenriched water at the site of cellulose synthesis; p ex is the proportion of exchangeable oxygen in carbohydrate substrates destined for cellulose synthesis; and ε bio is the biochemical enrichment (or fractionation) factor associated with the exchange of oxygen atoms between carbonyl group and medium water. The p x is expected to be unity for stem cellulose (Cernusak et al., 2005;Kahmen et al., 2011), but for leaf cellulose it is expected to be less than unity (Helliker and Ehleringer, 2002;Barbour, 2007;Kahmen et al., 2011;Song et al., 2014b; but see Liu et al., 2017), or even has been assumed to be zero . The p ex and ε bio have been successfully constrained using different approaches to be 0.42 ± 0.3 and 27 ± 3 , respectively (DeNiro and Epstein, 1981;Sternberg and DeNiro, 1983;Sternberg et al., 1986;Yakir and DeNiro, 1990;. These biochemical factors are a net result of complex biochemical processes in cellulose synthesis (Waterhouse et al., 2013) and seem to contain less climatic information, thus these factors are almost always treated as constants in mechanistic modeling. For instance, Helliker and Richter (2008) incorporated global tree-ring δ 18 O c data into mechanistic model and then inferred that boreal treecanopy temperature must be tremendously higher than ambient temperature to explain large observed leaf from their data. They also found that the inferred tree-canopy temperature throughout a large range of latitudes converges to ∼21 • C. This conclusion, however, was reached based on constant ε bio and p ex in their inverse modeling. Later, Sternberg and Ellsworth (2011) challenged this conclusion by showing that the ε bio is also temperature-dependent and increases with decreasing temperature. They heterotrophically generated cellulose from starch substrate at different temperatures and found that the ε bio could be as large as 31 when temperature is 5 • C. However, Zech et al. (2014a,b) indicated that the temperature-dependence of ε bio shown in that laboratory experiment might be an artifact from a variable p ex which, albeit not statistically significant (Sternberg, 2014), is temperature-dependent instead.
Here we use a global compilation of Sphagnum (peat mosses) δ 18 O c data to test the hypothesis that the ε bio is temperaturedependent. Sphagnum mosses are the keystone bryophyte species in peat-forming ecosystems, and these mosses are adapted to cool, waterlogged, and nutrient-poor conditions. Almost half of the boreal peatland areas are covered and built by Sphagnum mosses (Rydin and Jeglum, 2013). There have been a number of studies exploring the potential of stable isotopes in Sphagnum as proxies for peat-based paleoclimate reconstruction (e.g., Ménot-Combes et al., 2002;Daley et al., 2010;Loader et al., 2016;Granath et al., 2018;Xia et al., 2020). In general, Sphagnum mosses have many unique physiological characteristics with which their tissue δ 18 O c could closely track δ 18 O p with an offset dictated by ε bio . First, Sphagnum mosses are important peat-forming plants in bogs where moisture is predominantly derived from precipitation. Second, Sphagnum mosses have no stomata and no root systems. They maintain a simple water use strategy without the ability to control evapotranspiration. The water loss on Sphagnum surface is balanced either by new precipitation input or by capillary movement of water from below (Rydin and Jeglum, 2013;McCarter and Price, 2014). Therefore, Sphagnum mosses passively utilize water stored in surface unsaturated and oxic peat known as the acrotelm without the ability to access water at depth as vascular plants do. Third, Sphagnum mosses have the extraordinary capacity to store water and, more importantly, could protect the metabolic water from evaporation loss (Xia et al., 2020). Sphagnum have specialized large hyaline (water-filling) cells on their leaves and outer cortex of stems, and in some species these hyaline cells completely enclose photosynthetic cells (Loisel et al., 2009). Most of the water is, however, external capillary water held in pore spaces between leaves (Murray et al., 1989;Rydin and Jeglum, 2013). As a result of these multiple layers of water storage in Sphagnum mosses, the metabolic water in photosynthetic cells of Sphagnum living in local wet habitats could be isolated and weakly impacted by evaporative enrichment of 18 O (Xia et al., 2020). Fourth, Sphagnum mosses are opportunist photosynthesizers, growing and synthesizing new cellulose as long as conditions allow even in near-freezing temperature. This physiological trait allows us to explore the sensitivity of ε bio to air temperatures that are too low to allow the growth of vascular plants (Clymo and Hayward, 1982). In summary, the newly compiled dataset and analysis presented in this study offer a simple but effective test of the relationship between ε bio and temperature, while circumventing the complex physiology of vascular plants.

Sphagnum Oxygen Isotope Data
We compiled all available Sphagnum tissue δ 18 O measurements reported in literature, supplemented with some new measurements to fill regional gaps. The final dataset includes 785 individual δ 18 O measurements from 173 peatland sites (Figure 1). Key site information and references are summarized in Table 1, while raw data and metadata are in the Supplementary Dataset 1. Most of these sites are bogs from boreal and sub-arctic regions, but a few bog sites are located in the tropical and subtropical regions or the Southern Hemisphere. A few fen sites were also included in the dataset. The main Sphagnum species are the common bog species, S. magellanicum, S. fuscum, S. cuspidatum, and S. capillifolium. In some studies, the local microtopography information (such as hummock or hollow), or water table depth (WTD) for collected Sphagnum samples were available.
The exact materials of Sphagnum mosses used for isotope analysis differ among studies and include "capitula" (the cluster of new branch leaves at the tip of moss), "capitula cellulose, " "newest 1-cm stem cellulose, " "leaf, " "leaf cellulose, " "stem cellulose, " "whole plant cellulose, " and "5-cm whole plant cellulose." On the basis of our pilot experiments and literature data, we applied an offset correction to raw δ 18 O data that were measured from noncellulose ("capitula" and "leaf "), or stems ("newest 1-cm stem cellulose" and "stem cellulose"), and re-scaled these raw δ 18 O data to their equivalent capitula or leaf δ 18 O values with the involved uncertainty propagated in the following analysis. Correction is not required for δ 18 O data that were measured from "whole plant cellulose" or "5-cm whole plant cellulose" because our pilot experiment (n = 6) found that leaf tissues represent 90.3 ± 5.0% (error is the standard deviation, 1σ) dry mass of Sphagnum whole plant and have a similar cellulose yield (45.6 ± 4.7%) as stem tissues (48.8 ± 8.0%). The details about our cellulose extraction method can be found in Supplementary Text 1.

Precipitation Oxygen Isotope Data and Climate Data
We used the Online Isotopes in Precipitation Calculator (OIPC; Bowen, 2020) to obtain monthly δ 18 O p for every Sphagnum sampling location in the dataset. In the OIPC, δ 18 O p at a given location is modeled based on latitude and altitude plus a residual term that is interpolated using existing δ 18 O p observations (Bowen et al., 2005). To our knowledge, this model provides the best estimate for the site-specific δ 18 O p for our synthesis as there is no on-site precipitation δ 18 O measurement prior to Sphagnum sample collection except at three sites. The uncertainty of OIPC is depicted as the confidence interval for mean annual δ 18 O p prediction and is large in regions where δ 18 O p observations are scarce (Bowen and Revenaugh, 2003), but the confidence interval for monthly δ 18 O p prediction is not available in the OIPC. In addition, the OIPC is for predicting long-term mean value for monthly and annual δ 18 O p at a given location but does not account for their inter-annual variability. This inter-annual variability is an important uncertainty source in determining ε bio for Sphagnum as these samples were collected in different years and in some studies only their recent growth was harvested. To estimate this uncertainty, we used the archived outputs from three nudged isotope-enabled global circulation models (IsoGSM, LMDZ4, and GISS ModelE; accessible from https:// data.giss.nasa.gov/swing2/swing2_mirror/; Risi et al., 2012) to calculate the multi-model mean standard deviation of monthly δ 18 O p for each sample location. Therefore, we used the OIPC model that had been validated against observational δ 18 O p data to derive the site-specific δ 18 O p but used the isotope-enabled global circulation models nudged on reanalysis data to estimate the uncertainty associated with its inter-annual variability (at least for 31 years). We acknowledge that the accuracy of monthly δ 18 O p generated in the OIPC is unknown, but we used the prediction error for mean annual δ 18 O p as a reference.
We used the NASA MERRA-2 monthly 2-m air temperature (T2M in M2IMNXASM; GMAO, 2015b) data to infer the FIGURE 1 | Global map of peatland regions (green areas; Yu et al., 2010) and sites with Sphagnum δ 18 O c data presented in this study. Sites are color-coded to represent different regions: Alaska (blue), Canada and Contiguous US (orange), Europe (green), New Zealand (yellow), Patagonia (pink), Temperate East Asia (purple), Tropics and Subtropics (brown), and West Siberia (red).
Frontiers in Earth Science | www.frontiersin.org relevant growth temperature for each Sphagnum sample that was collected in different year and month. The MERRA-2 is an atmospheric reanalysis product for satellite era (beginning in 1980) and has a resolution of 0.5 • for latitude and 0.625 • for longitude. The raw temperature data were corrected for the altitude difference between each sampling site and the corresponding MERRA-2 grid assuming a lapse rate of −0.0065 • C/m. The MERRA-2 grid altitude is inferred from the constant model parameter surface geopotential height (PHIS in M2C0NXASM; GMAO, 2015a) as PHIS/g, where g is the standard gravity, 9.8 m/s 2 . Similarly, we also used the MERRA-2 monthly precipitation amount (PRECTOTCORR in M2TMNXFLX; GMAO, 2015c) data that were combined with OIPC δ 18 O p data to calculate the metrics of amount-weighted δ 18 O p that represent the isotopic composition of precipitation input for Sphagnum. However, the MERRA-2 products do not provide the uncertainty of reanalysis data other than the temporal variance of monthly mean.

Data Analysis
We first calculated the range and standard deviation of site-specific δ 18 O c data in one-time sampling campaign to characterize the within-site variability and evaluated if measured WTD or microtopographical locations would explain the within-site δ 18 O c variability. Then we calculated the apparent enrichment factor (ε app ) for each δ 18 O c datum as δ 18 O c minus its corresponding δ 18 O p , with the latter value depending on the type of materials and sampling time. We also similarly determined the relevant growth temperature for each sample from the MERRA-2 dataset. For δ 18 O c data derived from the recent growth of Sphagnum (capitula or the newest 1-cm stem), we used amountweighted δ 18 O p of the sampling month and the previous month (δ 18 O p2 ) to pair with temperature of the sampling month (T 1 ) for calculating the ε app . We assumed that the precipitation input in these recent 2 months was relevant to imprint the δ 18 O c signal in the tissue of recent growth, but growth was only restricted to the course of recent 1 month. To evaluate the sensitivity of calculated ε app to the choice of relevant time span for δ 18 O p and temperature, we also used δ 18 O p of the recent 1 month or 3 months (amount-weighted), together with the arithmetic average temperature of the recent two or 3 months as a pair. For δ 18 O c data derived from the whole plant of Sphagnum, we used the amount-weighted δ 18 O p of the entire recent 12 months (δ 18 O p12 ) to pair with the arithmetic average temperature of the most recent growing season defined as the months with MERRA-2 temperature above 0 • C during the recent 12-month period. We assumed that precipitation input during all these months either as rainfall or snowfall were relevant to imprint the δ 18 O c signal in whole plant, but growth itself was only restricted to the course of growing season. To minimize, but not necessarily eliminate the influence of evaporative enrichment that has caused the within-site variability in δ 18 O c and ε app , we used the sitespecific minimum δ 18 O c to calculate the site-specific minimum ε app , i.e., ε appmin , for each site as an approximation for ε bio . The precision of ε appmin was assessed using the Gaussian error propagation combining the analytical uncertainty in δ 18 O c , the uncertainty related to offset correction from raw δ 18 O data, and the climatological uncertainty in amount-weighted δ 18 O p values.
The accuracy of ε appmin was assessed using the OIPC prediction error for mean annual δ 18 O p . Pairing the calculated ε appmin and temperature data, we explored their relationship using a second order polynomial regression. Finally, we compared the result with previous study and discussed the implications for paleoclimate and plant physiology studies.

RESULTS
Our compiled Sphagnum tissue δ 18 O data ranged from 11.5 to 26.6 (relative to VSMOW). The uncertainty of δ 18 O analysis was represented by either the long-term precision of measurements (1σ) in different laboratories or the standard deviation (1σ) of replicate measurements, which was less than 0.6 and 2.0 , respectively. The difference in δ 18 O between extracted cellulose and untreated bulk tissue of the same Sphagnum material was small, and we determined the offset to be 0.69 ± 0.83 for capitula (error is the standard deviation, 1σ; Figure 2A). Kaislahti Tillman et al. (2010) reported that this offset is a negative value (Figure 2A), likely due to the fact that they measured Sphagnum materials from peat rather than from growing plants while untreated bulk tissues in peat probably had been contaminated by other chemical components in peat. The difference in δ 18 O between leaf cellulose (or capitula, which is formed by branch leaves) and stem cellulose of the same Sphagnum material was also small and well constrained in multiple studies as 0.62 ± 0.51 ( Figure 2B). This tissuespecific offset might be due to different water pools at the site of leaf and stem cellulose synthesis (Moschen et al., 2009). The capitula medium water might be subtly more enriched in 18 O than the matrix water below capitula where stem cellulose is synthesized. After offset correction, the median within-site ranges of δ 18 O c data were 1.0, 1.5,1.8, and 2.3 for sites with two, three, four, and at least five measurements, respectively ( Figure 3A). The within-site standard deviations (1σ) of δ 18 O c data were consistently between 0.7 and 0.8 , and are not influenced by the number of measurements per site ( Figure 3B). There is a pattern that more sites have a positive correlation between δ 18 O c and WTD ( Figure 3C) and that Sphagnum collected from hummock locations have higher δ 18 O c than from hollow locations ( Figure 3D), but few are statistically significant.
For all data from the recent growth of Sphagnum, there was a clear trend that the ε app increases with decreasing temperature, especially when T 1 is below 10 • C ( Figure 4A). The whole-plant data complemented the recent-growth data in the temperature range from 6 to 15 • C (Figure 4B), but their ε app did not show a clear trend with temperature. In both recent-growth and wholeplant data synthesis, outliers were mainly those data collected from fens and their ε app could be either higher, or lower than the rest of dataset. Unlike bogs that are rain-fed, fens are mainly fed by groundwater and stream water that that could complicate the source water δ 18 O for Sphagnum. In addition, Sphagnum inhabiting fens are not under desiccation stress and could develop different physiological strategy compared to those in bogs. These fen data were thus not included in the following analysis.
Merging both recent-growth and whole-plant data from bogs and plotting the site-specific ε appmin for sites with at least   Taylor (2008) and Brader (2013), plus a site from this study, in which samples were collected at microtopographical locations of hummock and hollow. Data from different sites are separated by horizontal lines. The "ns" denotes "not significant" difference at p < 0.05 level in Student's t-test between hummock and hollow δ 18 O c data for sites with at least two data points for both hummock and hollow. three data points versus corresponding growth temperature, the overall relationship that the ε appmin increases with decreasing temperature was slightly different from the regression models by Sternberg and Ellsworth (2011; Figure 4C). The second order polynomial regression from our data indeed corroborated the trend that the increase in ε appmin was steeper when temperature was lower, but the sensitivity was very weak for temperature range from 10 to 15 • C. Our additional analysis indicated that this pattern was not affected by pairing different time spans of precipitation input and temperature for recentgrowth dataset (Supplementary Figure 1). The ε appmin data points were considerably scattering relative to the trend line, but the deviations were acceptable considering the propagated uncertainty range. Notably, the accuracy of ε appmin was robust compared to their precision, except in a few cases, as shown by comparing the two types of error bars ( Figure 4C).

Within-Site Variability in Sphagnum Cellulose Oxygen Isotopes
The δ 18 O c variability for a specific site reflects differential evaporative enrichment of 18 O in metabolic leaf water for (C) The biplot between the site-specific ε appmin and growth temperature for bog sites with at least three measurements. The 95% prediction interval of OIPC-based mean annual δ 18 O p at each sampling location is shown as error bar with arrow caps. The precision (2σ) of ε appmin is denoted as error bar without caps. These data were fitted by a second order polynomial regression shown by the dash-dotted line. Also shown in these figures are the second order polynomial regressions fitted by data from heterotrophic generation of cellulose and data from submerged aquatic plants by Sternberg and Ellsworth (2011) as well as their extrapolations. Note that the full range of temperature (to 30 • C) in laboratory experiments by Sternberg and Ellsworth (2011) is shown to illustrate that the sensitivity of ε app to temperature is very weak at high temperature.
Sphagnum inhabiting different microtopographical locations in peatlands, as samples from the same site should receive the same precipitation input (Xia et al., 2020). This evaporation component-which could theoretically be expressed as leaf (1-p x p ex ) in Eq. (1)-has modified the δ 18 O p signal before the precipitation source water was incorporated into cellulose synthesis. For most sites, the δ 18 O c data were not measured from the same Sphagnum species, but recent studies suggested that the species effect of leaf anatomy or physiology (biotic factor) does not directly cause the within-site δ 18 O c variability (El Bilali and Patterson, 2012; Xia et al., 2020). Instead, different Sphagnum species favor different microhabitats (abiotic factor) such as drier hummocks, or wetter hollows that indirectly control the degree of evaporative enrichment of 18 O in metabolic leaf water and thus the δ 18 O c signal (Xia et al., 2020).
The result that most sites, if WTD data were available, displayed positive correlations between δ 18 O c and WTD suggests that Sphagnum samples associated with a deeper water table would register a higher degree of evaporative enrichment, although few of these correlations were statistically significant ( Figure 3C). Similar result was also found for data collected at hummock and hollow locations within a site (Figure 3D), despite that the actual WTD measurements and the size of these microtopographical features were unknown. The mechanism behind the connection between evaporative enrichment of 18 O in leaf water and WTD is more complex than a simple corollary that drier habitats with deeper water table cause a higher evaporation rate of leaf water and a higher enrichment of 18 O. In fact, the metabolic leaf water is technically difficult to collect as Sphagnum and other bryophytes store a large proportion of water in external capillary spaces (Dilks and Proctor, 1979). Through measuring Sphagnum water content and "bulk" leaf water δ 18 O, our previous study proposed that water stored in external capillary spaces of Sphagnum would evaporate first, thereby indirectly protecting internal metabolic leaf water from evaporative loss. As such, drier Sphagnum mosses would have a smaller proportion of external capillary water and as a consequence the metabolic leaf water is also subject to some degree of evaporation (Xia et al., 2020). That being said, the negative correlations between δ 18 O c and WTD also existed in the dataset (Figure 3C), although in at least a few cases these negative correlations were probably an artifact from a narrow range of WTD that did not represent a well-defined moisture gradient (Supplementary Figure 2). In addition, at a few sites Sphagnum inhabiting hollow could have higher δ 18 O c than hummock (Figure 3D). These opposite observations might be due to the decoupling between Sphagnum water content and WTD. Hollow-inhabiting Sphagnum have a higher collective surface area and these Sphagnum are more often susceptible to desiccation (Turetsky et al., 2008;Rydin and Jeglum, 2013) and could have a lower water content despite a shallower WTD, in particular after a sustained period of lacking precipitation. Studies that reported these data did not measure the Sphagnum water content nor had any description of their microhabitats, thus this specific hypothesis could not be directly evaluated.
Unlike vascular plants for which leaf is governed by leafair interaction through stomata and is highly dependent on relative humidity and transpiration rate (Zech et al., 2014a), there is no appropriate physiological model to mechanistically derive the leaf in Eq. (1) for Sphagnum mosses. Instead, the expression of leaf in Sphagnum should be considered in terms of the partitioning between internal and external leaf water rather than the atmospheric humidity. This water storage partitioning could explain why Sphagnum δ 18 O c data are less variable at site-specific scale ( Figure 3B) compared to their measured "bulk" leaf water δ 18 O that could be highly variable or highly enriched in 18 O (Ménot-Combes et al., 2002;Price et al., 2009;Loader et al., 2016;Xia et al., 2020). This "water buffer" mechanism might be unique for Sphagnum mosses adapted to waterlogged peatlands, but not for other bryophytes (Xia et al., 2020). For example, the moss species (Polytrichum strictum and Chorisodontium aciphyllum) from aerobic moss peatbanks in the Antarctic Peninsula showed very high and highly variable δ 18 O c values (Royles et al., 2013(Royles et al., , 2016, which are only possible with an unbuffered metabolic leaf water. Due to the presence of external water buffer for Sphagnum mosses, we propose that the evaporative enrichment of 18 O in their metabolic leaf water, despite being difficult to be directly measured or modeled, is a fairly small quantity compared to ε bio . From a process-based understanding, the evaporative enrichment should approximate to zero for Sphagnum inhabiting very wet habitats or having persistent external water layer, in which case the metabolic leaf water is close to being completely protected by external capillary water from evaporation. For individual sites where Sphagnum δ 18 O c data displayed a considerable variability, the site-specific minimum δ 18 O c value represents the case that the metabolic leaf water is least modified by evaporative enrichment. Therefore, we hereby justify that the corresponding ε app , i.e., ε appmin , serves as an alternative measure of ε bio for each of those sites. Future studies should aim to further refine the processes in how the isotopic signals of diverse Sphagnum water pools (e.g., Ooki et al., 2018) are modified evaporatively or transferred into newly formed cellulose in controlled laboratory experiments (Price et al., 2009;Brader et al., 2010).

Uncertainty in Constraining the Time Span of Precipitation Input and Sphagnum Growth
Defining the time span of growth and precipitation input for collected Sphagnum samples presents a main challenge in refining the relationship between ε app and temperature. The whole plant of Sphagnum likely integrated the growth for the time span of entire growing season, which is best defined as those unfrozen months (above 0 • C). The threshold value of 0 • C has been used to define bioclimatic variables such as the photosynthetically active radiation during the growing season (Gallego-Sala et al., 2010;Loisel et al., 2012;Charman et al., 2013). We thus used the T gs , the arithmetic average temperature of the months above 0 • C, as the growth temperature for the whole plant of Sphagnum. In fact, Sphagnum growth during the months of freezing temperature does occur if snow cover was thick and persistent enough to insulate Sphagnum from winter frost (Dorrepaal et al., 2004;Genet et al., 2013), as snow cover shading is not is not a major limiting factor for Sphagnum growth (Clymo and Hayward, 1982;Küttim et al., 2020). That being said, biomass production of Sphagnum during freezing months is much reduced. Considering this, using the T gs rather than the average temperature of every month of a year, could avoid underestimating the growth temperature for Sphagnum if the freezing season was very long or cold. However, the cumulative temperature-such that is defined as growing degree days above zero in bioclimatic studies (Charman et al., 2013)-might be a more accurate measure of the total warmth Sphagnum have experienced. Studies found that Sphagnum productivity at a given site is enhanced during the warmest season (Krebs et al., 2016;Küttim et al., 2020), although a few studies suggested that Sphagnum growth could be limited by summer desiccation stress and photoinhibition (Lindholm, 1990;Dorrepaal et al., 2004;Genet et al., 2013). Certainly, there is a tendency that T gs would underestimate the mass-weighted growth temperature for the whole plant of Sphagnum.
Although the freezing season is not considered relevant for Sphagnum growth, the accumulation of winter precipitation as snowpack during the freezing season could be an important component of source water for Sphagnum for the following snowmelt season. Indeed, direct water sampling of peatland water table showed δ 18 O values close to amount-weighted mean annual δ 18 O p (Ménot-Combes et al., 2002;Xia et al., 2020). For δ 18 O c data measured from the whole plant of Sphagnum, the amount-weighted δ 18 O p12 , which should be the most relevant representation for precipitation input, was used to derive the ε app to pair with the T gs in our analysis.
Unlike the whole plant, the capitula and the newest 1-cm stem only represent the recent growth of Sphagnum right before sample collection, and these tissues could register the recent input of δ 18 O p in cellulose. For example, Daley et al. (2010) collected Sphagnum capitula every two to three months and found that the δ 18 O c signals showed variations of more than 4 that were in general tracking the variations of δ 18 O p over seasonal scale at two of three sites in northwestern Europe, but at another site the variations were less than 2.5 , despite similar δ 18 O p seasonality at these locations according to the OIPC model. These observations imply that Sphagnum capitula regenerate and refresh their δ 18 O c signal at bimonthly and quarterly timescales, but the rate of δ 18 O c renewal for capitula appears to be variable. This isotopic renewal rate may depend on biotic factors such as Sphagnum productivity, which has been routinely determined by measuring Sphagnum stem height growth (Clymo, 1970). Field measurements indicate that stem height growth is enhanced with a higher photosynthetically active radiation and a longer growing season length (Loisel et al., 2012) but is also influenced by local ecohydrological factors (Gunnarsson, 2005;Loisel and Yu, 2013). The capitula turnover rate might be proportional to the stem increment rate as the capitula density is lower when stem growth is faster (Breeuwer et al., 2008;Loisel et al., 2012). Aldous (2002) used the ratio of number of branches per capitulum to number of branches per 1-cm of stem to normalize the turnover of capitula to stem increment. This method showed that capitula could replace themselves completely with every 1.65-cm stem increment for a case study of S. capillifolium from peat bogs in the northeastern United States (Aldous, 2002). Nevertheless, the exact time spans for the capitula replacement and the growth of newest 1-cm stem in the dataset are unknown. Other than the physical replacement, capitula δ 18 O c renewal rate may also depend on abiotic factors such as the residence time of metabolic leaf water δ 18 O, which in principle should be governed by both the size of acrotelm water reservoir and the rate of precipitation and evaporation. The role of residence time has been proposed to explain the smaller range of along-stem δ 18 O c variation compared to δ 18 O p variation in a case study of S. magellanicum collected from a southern Patagonian peat bog (Xia et al., 2020). The δ 18 O c renewal not only influences how many months of precipitation input are relevant for Sphagnum cellulose synthesis but also the average growth temperature, which are two important but potentially sitespecific metrics for our analysis on the temperature-dependence of ε app .
In our analysis of δ 18 O c data measured from the recent growth of Sphagnum, we used the δ 18 O p2 as the precipitation input to derive the ε app that was paired with the T 1 as the relevant growth temperature. This choice assumed that the sampling month temperature was the most relevant for capitula growth and the newest increment of stem, while metabolic leaf water not only derived from the precipitation of the sampling month but also partly inherited the precipitation input in the previous month. We notice that changing the definition of these time spans does not affect the overall trend between ε appmin and temperature, except in the case that the relationship between ε appmin and temperature became more like linear when the recent 3-month amount-weighted δ 18 O p was paired with the recent 3month average temperature (Supplementary Figure 1). Another possible complexity is that Sphagnum from the colder regions could have a slower capitulum turnover or stem increment rate (Loisel et al., 2012). If so, a longer time span of both δ 18 O p and temperature should be used specifically for these data. However, we argue that the site-specific local ecohydrological factors are equally important to climate factors to affect Sphagnum growth. For example, the growth rate of Sphagnum in cold arctic regions with very short summer could be remarkably high (Sonesson et al., 1980). That being said, when the recent 3-month amountweighted δ 18 O p and the recent 2-month average temperature were used instead for samples with T 1 lower than 6 • C, we found that the overall trend between ε appmin and temperature was not affected (Supplementary Figure 3).

The Effect of Temperature on Biochemical Enrichment Factor
The demonstrated relationship between Sphagnum ε appmin and growth temperature supports the previous finding that the ε bio is temperature-dependent in wheat seedling experiment by Sternberg and Ellsworth (2011). The relationship between ε bio and growth temperature in that experiment was established by assuming a known value of substrate δ 18 O and a constant of p ex , while the second order polynomial regression based on their experimental data was almost identical to that based on a compilation of aquatic plant δ 18 O c data (Sternberg and Ellsworth, 2011). From both their and our studies, the influence of temperature on ε bio is mainly manifested over the range of low temperature. Another previous study showing that the ε bio is insensitive to temperature  could be reconciled by the fact that the growth temperature for their plants was above 20 • C, at which the temperature sensitivity of ε bio is too weak to be convincingly detected.
Our results imply that the temperature-dependence of ε bio shown in heterotrophic generation of cellulose is not an artifact from a variable p ex as previously suggested (Zech et al., 2014a,b), because we used the site-specific minimum Sphagnum δ 18 O c data by which the term leaf was minimized and temperaturedriven variation in p ex , if any, was too small to introduce variation in δ 18 O c in Eq. (1). The leaf is usually a large quantity for vascular plants with stomata and increases substantially with decreasing relative humidity Zech et al., 2014a). Due to the large contrast between evaporated leaf water δ 18 O and unevaporated stem water δ 18 O in vascular plants, the p ex that could vary with plant types and in response to different growth conditions could play an important role in shaping the tree-ring δ 18 O c (Barbour and Farquhar, 2000;Ellsworth and Sternberg, 2014;Song et al., 2014a). Controversy has arisen from determining whether the p ex is dependent on temperature or not in laboratory experiments (Sternberg, 2014;Zech et al., 2014b), which in turn affects the determination of the corresponding ε bio . However, our study focusing on the field-collected Sphagnum mosses has avoided the complexity brought by the p ex . Finally, any variation in p ex could be further dampened by the factor of p x , which is less than unity for leaf tissue.
The second order polynomial regression from our dataset is slightly different from that by Sternberg and Ellsworth (2011) in that our regression line is flatter when temperature is above 10 • C and steeper when temperature is below 5 • C than their extrapolated line ( Figure 4C). We recognize that the growth temperature we used is essentially lower than the true temperature for moss growth that only occurs in warmer daytime. By comparison, both the temperatures for seedling growth and of aquarium medium water were controlled as experimental invariables and are more comparable to the daytime temperature for terrestrial plants in field conditions. Future studies should further investigate the discrepancy between mean temperature and photosynthetic temperature for naturally growing mosses and vascular plants.

Implications for Paleoclimate and Plant Physiology Studies
Plant δ 18 O c is a valuable paleoclimate proxy in a range of terrestrial climate archives, but studies employing plant δ 18 O c data should incorporate the influence of temperature on ε bio into the process-based proxy system model framework for interpreting proxy data (Dee et al., 2015;Keel et al., 2016;Okazaki and Yoshimura, 2019). The temperature effect might be small for Holocene-scale paleoclimate reconstruction, as temperature variability is in general less than 1 • C (Liu et al., 2014). However, this effect would be important for high-latitude and alpine regions as the ε bio exhibited the strongest sensitivity at low temperature and long-term temperature changes at these locations could be amplified due to various feedback mechanisms (Ballantyne et al., 2010;Feng et al., 2019). Furthermore, this temperature effect could be important in cold season, potentially resulting in varying biochemical enrichments in intra-annual tree-ring δ 18 O c records that could be used to reconstruct past climate seasonality (Schubert and Jahren, 2015;Zeng et al., 2016).
Plant δ 18 O c has also been used to understand the ecohydrological processes of plants such as tracing the water source for plants (Sargeant and Singer, 2016) or inversely modeling their growth temperature (Helliker and Richter, 2008;Song et al., 2011). Assuming a constant value of ε bio , Helliker and Richter (2008) showed that boreal tree-canopy temperature must be substantially higher than ambient temperature. Further analysis by Helliker et al. (2018) found that their previously inferred high tree-canopy temperature for boreal biome is essentially the result of a similarly warm daytime air temperature during growing season but tree-canopy itself could not reach such a high temperature. In fact, the temperature-dependence of ε bio is more important for subalpine and alpine regions where trees have a cool growing-season temperature (below 10 • C) and thus a considerably elevated ε bio . Indeed, existing tree-ring data from the Swiss Alps (Keel et al., 2016) and the Rocky Mountains in the western United States (Helliker et al., 2018) have shown that incorporating a temperature-dependent ε bio into mechanistic modeling would produce a better fit for observed δ 18 O c . Furthermore, in a study from the trans-Himalayan region (Chhumuja at an elevation of 4400 m above sea level) where growing season temperature is lower than 10 • C, a tree-ring δ 18 O c record was found to be best fitted with a relative humidity of 64%, despite the nearby station-observed values between 70% and 100% during the monsoon season (Brunello et al., 2019). This discrepancy could be alternatively explained by a higher ε bio at a low growing-season temperature.

SUMMARY
This study focused on Sphagnum mosses as a model plant and used an observational dataset of Sphagnum oxygen isotope measurements as a "natural experiment" to provide insights into the relationship between ε bio and growth temperature. After a careful consideration of Sphagnum physiology relevant for understanding the degree of 18 O enrichment in leaf water and the time span of precipitation input and plant growth, we demonstrated the temperature-dependence of ε bio for Sphagnum mosses comparable to what was determined experimentally in heterotrophic generation of cellulose. Although the relationship between photosynthetic temperature for field-growing mosses and mean air temperature remains to be refined and perhaps explains the discrepancy in two regression equations, both approaches suggest that the factor of temperature becomes increasingly important to increase the ε bio from the often assumed 27 when growth temperature is lower and thus should be considered in future paleoclimate and plant physiology studies wherever δ 18 O c data are measured, in particular for trees in high-elevation areas. Additionally, but not discussed here, our data compilation further constrained the within-site variability in Sphagnum δ 18 O c and implied that their δ 18 O c variations in peat core might reflect past changes in evaporative enrichment possibly linked to peatland surface patterning rather than δ 18 O p .

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online Supplementary Material.

AUTHOR CONTRIBUTIONS
ZX conceived the research, performed the laboratory and data analyses, and wrote the manuscript with input from ZY. ZY provided some samples and contributed to the discussion. Both authors contributed to the article and approved the submitted version.