Environmental and taxonomic controls of carbon and oxygen stable isotope composition in Sphagnum across broad climatic and geographic ranges

Rain-fed peatlands are dominated by peat mosses (Sphagnum sp.), which for their growth depend on nutrients, water and CO₂ uptake from the atmosphere. As the isotopic composition of carbon (^12,13C) and oxygen (^16,18O) of these Sphagnum mosses are affected by environmental conditions, Sphagnum tissue accumulated in peat constitutes a potential long-term archive that can be used for climate reconstruction. However, there is inadequate understanding of how isotope values are influenced by environmental conditions, which restricts their current use as environmental and palaeoenvironmental indicators. Here we tested (i) to what extent C and O isotopic variation in living tissue of Sphagnum is species-specific and associated with local hydrological gradients, climatic gradients (evapotranspiration, temperature, precipitation) and elevation; (ii) whether the C isotopic signature can be a proxy for net primary productivity (NPP) of Sphagnum; and (iii) to what extent Sphagnum tissue δ¹⁸O tracks the δ¹⁸O isotope signature of precipitation. In total, we analysed 337 samples from 93 sites across North America and Eurasia using two important peat-forming Sphagnum species (S. magellanicum, S. fuscum) common to the Holarctic realm. There were differences in δ¹³C values between species. For S. magellanicum δ¹³C decreased with increasing height above the water table (HWT, R²=17 %) and was positively correlated to productivity (R²=7 %). Together these two variables explained 46 % of the between-site variation in δ¹³C values. For S. fuscum, productivity was the only significant predictor of δ¹³C but had low explanatory power (total R²=6 %). For δ¹⁸O values, approximately 90 % of the variation was found between sites. Globally modelled annual δ¹⁸O values in precipitation explained 69 % of the between-site variation in tissue δ¹⁸O. S. magellanicum showed lower δ¹⁸O enrichment than S. fuscum (−0.83 ‰ lower). Elevation and climatic variables were weak predictors of tissue δ¹⁸O values after controlling for δ¹⁸O values of the precipitation. To summarize, our study provides evidence for (a) good predictability of tissue δ¹⁸O values from modelled annual δ¹⁸O values in precipitation, and (b) the possibility of relating tissue δ¹³C values to HWT and NPP, but this appears to be species-dependent. These results suggest that isotope composition can be used on a large scale for climatic reconstructions but that such models should be species-specific.

1 Introduction

Peatlands in temperate, boreal and arctic regions form large reservoirs of carbon, which are vulnerable to release under expected changes in global climate and land management (Rydin and Jeglum, 2013; Loisel et al., 2014). Because peat decomposes slowly and gradually accumulates, it preserves information on past peatland ecosystem dynamics and responses to allogenic and autogenic forcings. Palaeoenvironmental studies of peat may, therefore, help to anticipate the future responses of these globally important ecosystems to climate change (Loader et al., 2016). Past climate and local hydrology can be estimated using a variety of biotic and biogeochemical proxies, including the δ¹³C and δ¹⁸O values of organic material (e.g. van der Knaap, 2011; Royles et al., 2016; Lamentowicz et al., 2015). However, the environmental (e.g. climate) and biotic (e.g. species identity) controls of isotope differentiation in peatland-dwelling plants are still poorly understood and current assumptions regarding these controlling factors are yet to be tested on larger spatial scales.

Sphagnum mosses are the most dominant peat-forming plant group in acidic peatlands. The composition of stable isotopes of carbon and oxygen in their tissues is affected by different environmental conditions, operating through their impact on fractionation processes. When not submerged, carbon isotope signals in bulk tissues or components such as cellulose depend mainly on the concentration and isotopic composition of CO₂ in the chloroplast, which alters isotope discrimination during biochemical fixation of CO₂ (Farquhar et al., 1989; O'Leary, 1988). In mosses, the CO₂ concentration in the chloroplast, [CO₂]_c, is determined by temperature, light availability, CO₂ partial pressure and, most importantly, plant water status (Finsinger et al., 2013; van der Knaap et al., 2011; Ménot and Burns, 2001; Ménot-Combes et al., 2004; Royles et al., 2014; Skrzypek et al., 2007a; Kaislahti Tillman et al., 2013). When wet, external water films on leaf surfaces impede diffusion and [CO₂]_c is lowered (Rice and Giles, 1996; Rice, 2000; Williams and Flanagan, 1996); consequently, the proportion of fixed ¹³C increases due to internal drawdown of the preferred isotope ¹²C. When submerged, assimilation of respired or methane-derived CO₂ can alter [CO₂] and also the isotopic composition of C in Sphagnum (Raghoebarsing et al., 2005). Even when not submerged, respiratory carbon can be refixed by Sphagnum (Turetsky and Wieder, 1999; Limpens et al., 2008). Given that respired CO₂ is isotopically lighter than that in the atmosphere, it may also contribute to variation in tissue isotope values. Despite many detailed studies, there remains uncertainty about how the multiple controls on ¹³C isotope values combine to determine isotopic composition, and how universal the proposed mechanisms are on a global scale. This uncertainty currently restricts the utility of C isotope signals as palaeoclimatic or palaeoenvironmental indicators in peatlands (Loader et al., 2016).

Oxygen isotope values in moss tissues depend on the isotopic composition of the water sources, enrichment associated with evaporation from the moss surface and biochemical fractionation (Dawson et al., 2002). Once on the plant, ¹⁸O present in water equilibrates with that in atmospheric CO₂ prior to fixation as well as being incorporated directly during hydrolysis reactions, especially during the initial stages of carbon fixation (Gessler et al., 2014; Sternberg et al., 2006). Hence, variation in tissue oxygen isotopes reflect environmental conditions that control source water (rainfall, snowfall, groundwater) as well as fractionation caused by evaporation prior to fixation, which is controlled by micrometeorological conditions (mainly temperature, relative humidity and incident energy) (Daley et al., 2010; Moschen et al., 2009; Royles et al., 2013; Kaislahti Tillman et al., 2010). Oxygen isotope composition has, therefore, been used to reconstruct climatic conditions and to infer the dominant water source in peatlands (Aravena and Warner, 1992; Ellis and Rochefort, 2006; van der Knaap et al., 2011). Ongoing measurements of oxygen isotopes in precipitation across the globe (Bowen, 2010; IAEA/WMO, 2015) have generated models that predict spatial patterns in oxygen isotope composition of precipitation based on temperature, elevation, atmospheric residence time and circulation patterns (e.g. Bowen, 2010). Once isotopic composition of the source water is accounted for, variation in moss tissue isotopic values should be largely determined by fractionation that accompanies evaporation from the surface of plants. How well oxygen isotopes in Sphagnum tissues reflect atmospheric water or plant surface water depends on local weather conditions such as precipitation, air temperature and humidity. For example, Bilali et al. (2013) suggest that oxygen isotopes in Sphagnum mosses from maritime bogs will track variation in precipitation patterns whereas isotopic values in continental habitats will be more dependent on summer temperature, as temperature and humidity are more variable in those regions. On local scales, oxygen isotope values also vary as a function of temperature and humidity. Aravena and Warner (1992) found differences that correspond with changes in microtopography. Elevated microsites (hummocks) were enriched in ¹⁸O, which they ascribed to higher evaporation compared to that of neighbouring wet depressions (hollows). However, as with ¹³C, there remains uncertainty in how ¹⁸O signatures relate to environmental factors and species identity and to what extent global ¹⁸O patterns in precipitation dominate over local processes.

Stable isotopes can also serve as indicators of net primary productivity (NPP) (Rice and Giles, 1996; Williams and Flanagan, 1996; Rice, 2000). However, few studies have explored these relationships in the field. In a multispecies comparison of peat mosses, Rice (2000) found that plants with higher relative growth rates had lower discrimination against ¹³C and therefore were more enriched in ¹³C. This was attributed to the local environment, with fast-growing plants of wetter microhabitats having thicker water films that inhibit CO₂ diffusion into the plant and to species-specific differences in maximum rates of photosynthesis. Both factors would reduce internal [CO₂] and thereby lower discrimination. In line with this, a warming experiment by Deane-Coe et al. (2015) reported a positive relationship between moss net NPP and δ¹³C values for tundra mosses (Dicranum, Pleurozium, Sphagnum). Clearly, carbon isotope values show promise as indicators of peat moss contemporary growth and potentially as an NPP proxy in palaeoecological studies. This could be particularly valuable to differentiate productivity and decomposition controls in long-term carbon accumulation studies. To date, we are not aware of attempts to explore the robustness of these relationships across large spatial scales.

Together, tissue carbon and oxygen isotope compositions are controlled both by environmental factors at micro- and macro-scales, and by species-specific differences that relate to water balance and carbon dynamics in peat mosses. Palaeoecological studies rely on such environment–isotope relationships for environmental reconstructions (Ellis and Rochefort, 2006; van der Knaap et al., 2011). The underlying mechanisms are, however, rarely fully explored using known environmental gradients (but see Ménot and Burns, 2001 for an example) or only tested across narrow bands of environmental variation, often with sets of correlated environmental factors (Loader et al., 2016). Moreover, interactions with biotic factors such as species identity have received little attention despite the large variations in Sphagnum species dominance commonly observed down peat cores (e.g. Ménot and Burns, 2001). Here we aim to provide a robust, cross-scale evaluation of how environmental factors and species identity influence the C and O isotope compositions of Sphagnum using two common and widely distributed peat-forming species (S. magellanicum and S. fuscum) that are primarily rain-fed. To achieve this, we performed an unprecedented large sampling campaign across the Holarctic realm.

Specifically, we (i) investigated relationships between C and O isotope values and factors known to influence plant water availability (height above the water table – HWT, temperature, evaporation and precipitation) and CO₂ partial pressure (elevation), and tested whether their effects were modified by species identity; (ii) tested the prediction that Sphagnum tissue δ¹³C values are associated with NPP; and (iii) tested whether tissue δ¹⁸O in rain-fed Sphagna is predicted by the δ¹⁸O isotope signature in precipitation but modified by negative relationships with precipitation and positive ones with temperature/evaporation. Across these objectives we examined how C and O isotope values varied with scale (within-peatland versus between-peatlands) and to what extent HWT and NPP could explain variation within and between peatlands.

2 Materials and methods

2.1 Study species and collection sites

Our study focused on two common peat-forming Sphagnum species, S. fuscum (Schimp.) H. Klinggr. (circumpolar distribution) and S. magellanicum Brid. (cosmopolitan distribution). In general, these species are confined to primarily rain-fed peatlands (bogs) and described as hummock (S. fuscum) and lawn (S. magellanicum) species. However, S. magellanicum is a species with a very broad niche and found in a range of habitats with varying degrees of groundwater influence (Flatberg, 2013). These species are easy to identify, but recent research has shown that the dark European morph of S. fuscum is conspecific to the North American S. beothuk (Kyrkjeeide et al., 2015), and S. magellanicum has been shown to consist of two genetically diverged morphotypes (Kyrkjeeide et al., 2016) that recently were separated at the species level (Hassel et al., 2018). Unpublished genetic data suggest that samples collected in our study consist of both S. magellanicum morphs (approximately 50/50) and possibly one or two samples of S. beothuk (Narjes Yousefi, personal communication, 2018). Hence, we here treat our species as aggregates (i.e. species collectiva), S. fuscum coll. and S. magellanicum coll.

The two species were sampled across the Holarctic region at a total of 93 sites (Fig. 1; Table S1 in the Supplement) at the end of the growing season. To make comparisons between species and between sites possible, we focused on habitats where both species can be found and have low influence of surrounding groundwater. Thus, we only sampled bogs (including a few poor fens with ombrotrophic character) and open (no tree canopy) habitats. Sampling was conducted mainly during 2013, but a few sites were sampled at a similar time of year in 2014. At each site two patches (minimum 10 m apart) for each species were sampled (except for 11 sites that contained only one patch for one species). At each sampling patch we recorded moss growth, HWT and GPS coordinates, and collected a moss sample (78 cm² and 5 cm deep) at the end of the growing season (September to November depending on location and generally coincided with when there was a risk of the first snowfall to occur). Moss samples were dried (24 h at 60–65 ^∘C) within 72 h or immediately frozen and later thawed and dried. The apical part (the capitula, top 1 cm) of the dried plant shoots was used for isotope analysis, while the stem section was used for bulk density estimation to calculate moss NPP.

Figure 1 Map illustrating sample sites for the investigated species. At some sites only one of the two Sphagnum species was sampled, indicated by red triangles or black circles. Otherwise sites contained both species (blue crosses). The map is centred on the North Pole and has an orthographic projection. Geographical ranges: latitude 41.6–69.1 N, elevation 2 – 1829 m a.s.l. See Table S1 for details.

2.2 Isotope determination

Ten capitula from each patch were selected and finely chopped with a single-edge razor by hand and mixed. Capitula were chosen as they reflect the most recently fixed organic matter and should relate better to recent growing season conditions. In Sphagnum, δ¹³C from the capitulum is similar to that of branches within the top 15 cm of plants but is approximately 1–2 ‰ less negative than stems (Loader et al., 2007). For δ¹⁸O, the offset between branches and stems is around 1 ‰ (Moschen et al., 2009). Standard deviations of repeated samples were 0.6 ‰ and 0.7 ‰ for δ¹³C and δ¹⁸O, respectively. Approximately 0.5 mg of dry sample was packed in tin cups for δ¹³C analyses, and ∼0.2 mg in silver cups for δ¹⁸O analyses. Samples were analysed at Union College (Schenectady, NY, USA) using a Thermo Delta Advantage mass spectrometer in continuous flow mode connected (via a ConFlo IV) to a Costech Elemental Analyzer for δ¹³C analysis or a Thermo TC/EA for δ¹⁸O analyses. Isotope values are presented as $\mathrm{1000}\times (R_{\mathrm{sample}}/R_{\mathrm{standard}}-\mathrm{1})$, where R_sample and R_standard are the ratios of heavy to light isotopes (e.g. ¹³C ∕ ¹²C) and are referenced to VPDB and VSMOW for C and O, respectively. Carbon isotope data were corrected using sucrose (IAEA-CH-6, −10.449 ‰), acetanilide (in house, −37.07 ‰) and caffeine (IAEA-600, −27.771 ‰). Oxygen isotope data were corrected using sucrose (IAEA-CH-6, 36.4 ‰), cellulose (IAEA-C3, 31.9 ‰) and caffeine (IAEA-600, −3.5 ‰) with values from Hunsinger et al. (2010). Oxygen isotope standardization was further checked with the whole wood standards USGS54 and USGS56. The combined instrument uncertainty for δ¹³C (VPDB) is <0.1 ‰ based on the in-house acetanilide standard and <0.5 ‰ for δ¹⁸O (VSMOW) based on the cellulose standard (IAEA-C3).

We performed isotope analyses on whole-plant tissue rather than on cellulose extracts. In living Sphagnum samples, there is a strong linear relationship between the isotopic composition of these two components for both δ¹³C (R² values 0.89–0.96; Kaislahti Tillman et al., 2010; Ménot and Burns, 2001; Skrzypek et al., 2007b) and for δ¹⁸O (R² values 0.53–0.69; Kaislahti Tillman et al., 2010; Jones et al., 2014). Focussing on whole-plant tissue allowed us to analyse a higher number of samples for this study, allowing larger numbers of sites and more replication.

Table 1 Sample sizes, standard variation and overall partitioning of measured variation for each species and response (δ¹³C and δ¹⁸O). N_site is the number of sites and N_obs the total sample size. Standard deviation of the responses is given for within and between sites, together with the proportion of total variance measured between sites and within sites.

2.3 Environmental variables

The modelled δ¹⁸O signal in meteoric water (precipitation) (Bowen and Wilkinson, 2002) was obtained from http://www.waterisotopes.org (last access: 2 October 2017) as annual and monthly isotope ratio estimates at 10 arcmin resolution. These global estimates have shown to be highly accurate (R²=0.76 for mean annual δ¹⁸O in precipitation) and are based on absolute latitude and elevation and account for regional effects on atmospheric circulation patterns (for details see Bowen, 2010, 2017; IAEA/WMO, 2015). To test which temporal period of δ¹⁸O values in precipitation showed the highest correlation with tissue δ¹⁸O values, we calculated annual (January–December), growing season (May–October) and winter–spring (January–April) mean isotope ratio. We calculated both unweighted and weighted means against precipitation for each month. Monthly precipitation (PRECTOTCORR), land evapotranspiration (EVLAND) and surface air temperature (TLML) for each site and year of sampling (2013 or 2014) were retrieved from the NASA GESDISC data archive, land surface and flux diagnostics products (M2T1NXLND, M2TMNXFLX; resolution longitude 0.667^∘, latitude 0.5^∘; Global Modeling and Assimilation Office, 2015a, b). Total precipitation and evapotranspiration (ET), and mean temperature, from April to October were used as predictors in the statistical models. As ET can be compensated for by precipitation, we used the ET/P quotient as a predictor for the effect of water loss. A high value (>1) indicates a net loss of water to the atmosphere. Site altitude was retrieved from a global database using the R package elevatr (ver 0.1-2, Hollister and Shah, 2017).

The distance from the moss surface to the water table (HWT) was measured using water wells (commonly a PVC pipe, 2–5 cm in diameter and slotted or perforated along the sides) with a “plumper” (a cylinder on a string that makes a “plump” sound when it hits the water surface) or a “bubbler” (a narrow tube that makes bubbles when it hits the water surface while the user blows into it). HWT was measured in the spring and in the autumn and there was a strong correlation between the two time points (r=0.74). As growth mainly occurs in late summer to autumn in temperate and boreal regions, we used HWT at the end of the season as the proxy of relative HWT between sites.

2.4 Moss growth

Moss growth (or productivity, NPP) was measured with a modified version of the cranked wire method (see Clymo, 1970; Rydin and Jeglum, 2013 for details), with bristles from a paint brush spirally attached to a wire. These “brush wires” were inserted in the moss layer with the end of the wire protruding above the surface. Height increment (i.e. vertical growth) was measured over the growing season as the change in distance (to nearest millimetres) between the moss layer and the top of the wire. A minimum of three wires were inserted within a 1×1 m uniform area (same microhabitat, vegetation and general structure). To determine moss bulk density (kg m⁻³), we dried (24 h at 60–65 ^∘C) the top 30 mm of the stems (area 78 cm²) in our collected core (see Sect. 2.1). Biomass growth on an area basis (g m⁻² yr⁻¹) was calculated as height increment × bulk density.

2.5 Statistical analyses

To test and quantify the influence of environmental variables and species identity on isotope composition, we used linear mixed models in R (R core team, 2016), employing the R package lme4 ver 1.1-12 (Bates et al., 2015). Site dependence (i.e. multiple samples from the same site) was accounted for by adding site as a random factor. For tissue δ¹³C, we first fitted two separate models to test the independent effects of HWT, NPP and species identity (S. fuscum and S. magellanicum), and whether the HWT or NPP effect varied between species by fitting a species interaction term. To test the explanatory power of environmental variables (ET/P, precipitation, temperature, elevation) we first constructed a base model with HWT and NPP, as they were identified as the main predictors in literature. For simplicity we removed negligible interactions from this model. Each environmental variable and its interaction with species was then tested against the base model. For tissue δ¹⁸O, we first explored which temporal period of modelled δ¹⁸O_precip (annual, growing season, winter–spring) had the highest explanatory power and whether the relationship varied between species. The identified best model was then used as the base model to separately test each environmental variable (HWT, ET/P, precipitation, temperature) and its interaction with the species.

The proportion of variance explained by the predictors was calculated at the site level (Gelman and Hill, 2007) or as marginal R² (Nakagawa and Schielzeth, 2013; R package piecewiseSEM ver 1.1-4; Lefcheck, 2016). Although our study focused on explained variance by predictors, we also performed statistical tests of predictors and their interactions using type-2 (main effects tested after all the others in the model but without the interaction term) F tests, applying Kenward–Roger adjustments to the degrees of freedom, as implemented in the car package (ver. 2.1-3, Fox and Weisberg, 2011). Standard model checking was performed (e.g. residual analyses and distribution of random effects) to ensure compliance with model assumptions. Covariances between predictors were small (r<0.15) or moderate (r=0.40–0.50 between ET/P, precipitation and temperature) and this multicollinearity had minor impact on model estimates.

Table 2 Results from linear mixed models for δ¹³C values. Statistical tests are based on type-2 F test using Kenward–Roger-adjusted degrees of freedom. The second model only included S. magellanicum. Elevation (m a.s.l.) and the three climatic variables (growing season sums and means: ET/P, P in millimetres, temperature in degrees Celcius) were tested one by one in the model including HWT (height above the water table in centimetres), species and NPP (g m⁻² yr⁻¹). For simplicity, the negligible HWT×NPP term was dropped from this model (P=0.36). Estimated effects (± SEs) are only given for the main effects if interactions were considered negligible. These effects are slopes for continuous variables (all variables except species) and for species (categorical) the difference between S. magellanicum and S. fuscum (i.e. S. fuscum being the reference level). In the presence of an interaction between HWT and species, the species effect was estimated at mean HWT. $R_{\mathrm{site}}^{\mathrm{2}}$ is explained between-site variance; $R_{\mathrm{marginal}}^{\mathrm{2}}$ is explained total variance.

^* The effect of S. magellanicum compared to S. fuscum at HWT 28 cm.

3 Results

Data collection from a geographically broad area resulted in large variation of isotope values and explanatory variables (Table S2). Due to uncertainty in height increment measurement we recorded a few negative values resulting in negative NPP. These values were kept in the analyses. The means and standard deviations (in parenthesis) of height increment (HI, millimetres) were 14.3 (10.1) for S. fuscum and 19.5 (14.1) for S. magellanicum, and the means and standard deviations of bulk density (BD, kg m⁻³) were 17.8 (9.9) for S. fuscum and 10.2 (7.6) for S. magellanicum.

Table 3 Results from linear mixed models for δ¹⁸O values. Statistical tests are based on type-2 F test using Kenward–Roger-adjusted degrees of freedom. Three time periods for modelled δ¹⁸O values (‰) in precipitation were tested individually: annual mean, growing season (April–September) and spring (January–April). The three climatic variables (growing season sums and mean: ET/P, P (mm), temp [^∘C]) were tested one by one in a model including HWT (cm) and mean annual δ¹⁸O values). Estimated effects (± SEs) are only given for the main effects if interactions were considered negligible. These effects are slopes for continuous variables (all variables except species) and for species (categorical) the difference between S. magellanicum and S. fuscum (i.e. S. fuscum being the reference level). $R_{\mathrm{site}}^{\mathrm{2}}$ is explained between-site variance; $R_{\mathrm{marginal}}^{\mathrm{2}}$ is explained total variance.

3.1 δ¹³C signal

Variation in Sphagnum tissue δ¹³C values was marginally greater within sites than between sites (Table 1). HWT predicted the δ¹³C values, but the relationship differed between the two species (Table 2, Fig. 2). Although δ¹³C values decreased with increasing HWT for both species, the slope was less steep for S. fuscum and this species had slightly higher δ¹³C values overall. In separate models for the two species, HWT for S. fuscum had near-zero explanatory power, while for S. magellanicum HWT explained 33 % of the between-site variation and 17 % of the total variance (i.e. marginal R²).

Figure 2 Relationship between height above the water table (HWT, measured at the end of the growing season) and δ¹³C values in two Sphagnum species sampled across the Holarctic realm. Lines show the pooled regression line in a mixed-effect model. Shaded areas indicate approximate 95 % confidence intervals (2 × SE) of the regression coefficients that do not include the uncertainty of the random effects. Equations: S. fuscum, δ¹³C = $-\mathrm{27.56}-\mathrm{0.021}\times \mathrm{HWT}$; S. magellanicum, δ¹³C = $-\mathrm{27.74}-\mathrm{0.045}\times \mathrm{HWT}$. N_site=83, N_total=311. See also Table 1.

Measured δ¹³C values were related to moss NPP and δ¹³C values increased by 0.0023 ‰ (SE: 0.00048) for each milligram of biomass produced per square metres. NPP explained 11 % of the between-site variation in δ¹³C and 7 % of the total variation. HWT and NPP explained 48 % of the between-site variation of δ¹³C in S. magellanicum and 24 % of the total variation. Corresponding values for S. fuscum were 6 % and 7 %. Of the additional environmental variables tested, we found weak evidence that ET/P and temperature were positively correlated with δ¹³C but only for S. magellanicum (Table 2).

3.2 δ¹⁸O signal

Sphagnum tissue δ¹⁸O values varied more between sites than within sites, and at similar magnitudes and proportions for both species (Table 1). Tissue δ¹⁸O values were predicted by the spatially explicit estimates of δ¹⁸O value isotope signature in precipitation (Fig. 3, Table 3). Annual mean δ¹⁸O_precip explained 69 % of the variation in δ¹⁸O_tissue between sites. This was similar to mean winter–spring (January–April) δ¹⁸O_precip values (75 % explained) but higher than the growing season (April–September) δ¹⁸O_precip (58 %). Using precipitation-weighted δ¹⁸O_precip values resulted in lower percentages of explained variance for all three time periods ($R_{\mathrm{site}}^{\mathrm{2}}$: annual 52 %, January–April 65 %, April–September 52 %). S. magellanicum had consistently lower δ¹⁸O values than S. fuscum (−0.83 ‰), but both species had a similar relationship between tissue δ¹⁸O and δ¹⁸O_precip (Fig. 3, Table 3).

Figure 3 Association between modelled annual mean δ¹⁸O_precip values and δ¹⁸O values in two Sphagnum species. Data show site means for each species and error bars represent the standard deviation (at some sites only one sample of a species was taken). Regression lines with different intercepts (P<0.001, Table 2) illustrate the relationship between modelled δ¹⁸O_precip and Sphagnum δ¹⁸O. Equations: S. fuscum, $\mathrm{26.36}+\mathrm{0.43}\times {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{precip}}$ (n=1–2, N_site=80); S. magellanicum, $\mathrm{25.53}+\mathrm{0.43}\times {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{precip}}$ (n=1–2, N_site=83). Shaded areas indicate approximate 95 % confidence intervals (2×SE) of the regression coefficients that do not include the uncertainty of the random effects.

HWT at the end of the growing season was, on average, 11 cm lower in S. magellanicum patches (wetter habitat) compared to S. fuscum (HWT = 33 cm) patches ($F_{\mathrm{1},\mathrm{224}}=\mathrm{131.9}$, P<0.0001). However, we found only very weak support for the hypothesis that HWT predicts tissue δ¹⁸O values, as HWT explained <1 % of the δ¹⁸O variation (Table 2). There was negligible influence of the additional environmental variables on δ¹⁸O values (Table 2). ET/P was associated with higher δ¹⁸O values in S. magellanicum and lower in S. fuscum (but not different from the zero effect), while increasing temperature was weakly associated with overall lower δ¹⁸O values.

4 Discussion

4.1 Stable carbon isotope discrimination in Sphagnum

Our data were consistent with the hypothesis that moss growing closer to the water table (low HWT) has reduced carbon isotope fractionation, leading to greater fixation of ¹³CO₂ and more ¹³C-enriched tissue (Rice and Giles, 1996; Williams and Flanagan, 1996). Given that the water table position was measured in different places at different times and all are one-time measurements, this result is remarkably robust. For example refixation of ¹²C-enriched substrate-derived CO₂ in living Sphagna (Turetsky and Wieder, 1999; Raghoebarsing et al., 2005) can potentially contribute to within-site variation in δ¹³C as it potentially affects both the ambient concentration of CO₂ as well as its isotopic composition. Interestingly, the strength of the δ¹³C – HWT relationship differed in the two species, with S. magellanicum exhibiting a greater reduction in δ¹³C in response to drier conditions (high HWT) than S. fuscum. The weaker effect of HWT on δ¹³C values in S. fuscum is likely a consequence of limited fluctuation in tissue water content, as this species is well known to store abundant water within capillary spaces and resist drying (Rydin, 1985), thus maintaining the waterfilm that results in reduced fractionation. Loader et al. (2016) reported a similar slope estimate for S. magellanicum in a single peatland, and several studies have confirmed the effects of contrasting microtopography (i.e. hummock – hollow differences) using multi-species comparisons (Price et al., 1997; Loisel et al., 2009; Markel et al., 2010). As such, our results suggest that species-specific differences in carbon isotope discrimination in Sphagnum are related to water retention capacity and, consequently, become more apparent under drier conditions. This supports the results of previous, smaller-scale studies (Rice, 2000).

The influence of species identity on the relationship between δ¹³C values and water table position has important implications for palaeoenvironmental reconstructions based on δ¹³C values. The relationship between δ¹³C and HWT has been used in palaeoecological reconstructions of surface wetness (e.g. Loisel et al., 2009). In our data set the strength of the relationship was weaker than previously reported. For instance, Loader et al. (2016) reported R²=54 % for S. magellanicum in a single site. Given the characteristics of our data (large-scale, circumpolar), the explanatory power ($R_{\mathrm{marginal}}^{\mathrm{2}}=\mathrm{17}$ %) can be considered acceptable and comparable to other proxies such as testate amoebae (16 % in Loader et al., 2016; Sullivan and Booth, 2011). Our results imply that isotopic signals of peatland wetness in hummock-dwelling species (such as S. fuscum) may be weaker or absent compared to lawn species. It is therefore important that the same species or species type (e.g. lawn species as they likely have a broad HWT niche) is used if δ¹³C values are employed as a proxy to infer changes in HWT.

We also identified evidence that evapotranspiration (ET) and NPP modify δ¹³C values, although the effect of ET was weak and restricted to S. magellanicum. We expected a stronger relationship, as ET and temperature control δ¹³C by increasing water loss at the moss surface and reducing the diffusive resistance (i.e. reducing CO₂ limitation), enabling increased discrimination against ¹³C (Williams and Flanagan, 1996). NPP only explained a small proportion of the variation in δ¹³C values, but the relationship was apparent across species. Several studies have proposed the use of δ¹³C values to infer Sphagnum productivity (e.g. Rice and Giles, 1996; Rice, 2000; Munir et al., 2017) and our study is the first to test this at the pan-Holarctic scale. Deane-Coe et al. (2015) investigated δ¹³C values across moss species (including Sphagnum) and years at one site and found a weak relationship between productivity and δ¹³C values (R²=0.10 and 0.31). Similarly, Rice (2000) reported that the relative growth rate explained about 25 % of the variation in ¹³C discrimination. We did not find as strong a relationship (R²<0.12), but our study was geographically broader and less controlled; consequently, our results were likely influenced by more complex interactions among environmental factors that affect Sphagnum growth across our sites. Nevertheless, our results indicate independent effects of evaporation and productivity on δ¹³C values. The lack of a strong NPP pattern somewhat limits the ability to infer productivity of Sphagnum in palaeoecological studies.

4.2 Global patterns of δ¹⁸O values in Sphagnum

Modelled δ¹⁸O values in precipitation (Bowen, 2010) explained much of the variation in δ¹⁸O_tissue values between sites (R²=68 % for annual mean δ¹⁸O_precip). The percent variance explained was even higher if the spring period for modelled ¹⁸O_precip was used but was lower for the growing season average. This result does not necessarily mean that spring season water was utilized by the plants during the growing season. Between-site variation in ¹⁸O_precip values are much larger in the winter (Fig. S1), more effectively discriminating maritime and continental regions (Bowen, 2010). The better fit may simply be an effect of a more distinct separation of ¹⁸O_precip in the winter data. Although the δ¹⁸O_tissue−¹⁸O_precip relationship presented here is robust, a few δ¹⁸O values are less well predicted by the regression model and they originate from Northwest Territories (Canada) and West Siberia (Russia). Likely, this suggests that the ¹⁸O_precip model is less accurate in these interior regions with fewer precipitation collection stations.

In contrast, the data did not support a negative correlation between precipitation amount and δ¹⁸O_tissue values and δ¹⁸O_tissue values were only weakly affected by predictors associated with water loss (ET/P and/or temperature) and species identity. The indication of ¹⁸O enrichment in S. magellanicum due to ET/P was expected as the lighter isotope ¹⁶O needs less energy to vaporize. However, the opposite trend was suggested for S. fuscum, and surprisingly, higher surface temperatures decreased ¹⁸O enrichment. Hence we conclude that climatic variables associated with water loss were weak predictors after controlling for δ¹⁸O_precip values. This result may not be too unexpected as laboratory experiments have so far failed to relate ¹⁸O enrichment in Sphagnum to differences in evaporation rates (Brader et al., 2010).

There have been few regional studies on moss δ¹⁸O_tissue values that span gradients of δ¹⁸O_precip values (Royles et al., 2016; Skrzypek et al., 2010) and most interpretations of moss δ¹⁸O_tissue – climate relationships come from peat core studies (e.g. van der Knaap et al., 2011). In Antarctic non-Sphagnum peat banks variation in δ¹⁸O_cellulose values tracked δ¹⁸O values in moss water across a latitudinal gradient (61–65^∘ S) despite a lack of difference in δ¹⁸O_precip. This result led Royles et al. (2016) to suggest that moss water and tissue δ¹⁸O values are better temporal integrators of source water than point rainfall measurements. The authors interpreted site-to-site differences as relating to differential evaporative enrichment and other physio-chemical factors that affect ¹⁸O exchange, fixation and biochemical synthesis. Skrzypek et al. (2010) explored variation in Sphagnum δ¹⁸O_tissue values across a regional altitudinal gradient and found no consistent trend or significant relationship linking δ¹⁸O_tissue values to altitude, where δ¹⁸O in source water is expected to differ. Although fractionation in source water caused by adiabatic cooling with altitude should lead to altitudinal effects, differences in precipitation amount can confound this pattern (Gat et al., 2000). Unfortunately, there are limited regional studies that have tested the effects of variation in source water on δ¹⁸O_tissue values. The present study provides a much greater range of geographical and environmental variation and shows strong support for Sphagnum strongly tracking source water.

Interestingly, the relationship between δ¹⁸O_tissue and δ¹⁸O_precip values detected here is very similar to that proposed some time ago by Epstein et al. (1977); ${\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{cellulose}}=\mathrm{27.33}+\mathrm{0.33}\times {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{precip}}$ (note that Jones et al., 2014, show high correspondence between δ¹⁸O_cellulose and δ¹⁸O_tissue values). However, our data suggest a slightly steeper slope and lower intercept, particularly for S. magellanicum. The species effect on δ¹⁸O suggests a difference in the degree of evaporation from the plant surface prior to the uptake of water. The lower δ¹⁸O values for S. magellanicum compared to S. fuscum (−0.83 ‰) is comparable to the results from bogs in Canada for the same species (−2.2 ‰, Aravena and Warner, 1992) and between a hollow and a hummock species in the Netherlands (−2 ‰, Brenninkmeijer et al., 1982). This suggests that the absorbed water in S. magellanicum was subject to less evaporation. In Sphagnum plants, surface water is largely affected by capillarity, water storage and reducing conductance with compact morphology. Plant traits that enhance these functions are more pronounced in species and individuals found at high HWT as these characteristics maintain high tissue water content (Hayward and Clymo, 1982; Laing et al., 2014; Waddington et al., 2015). Consequently, during droughts, Sphagnum species growing close to the water table will dry out quickly as the evaporative demand cannot be balanced, and simultaneously photosynthesis is shut down. Sphagnum species higher above the water table wick water from below and store water effectively, thereby remaining photosynthetically active while water is lost due to evaporation. This mechanism would result in ¹⁸O enrichment being higher above the water table (Brenninkmeijer et al., 1982; Aravena and Warner, 1992) and explains the positive relationship between HWT and δ¹⁸O in S. magellanicum reported by Loader et al. (2016) along a 10 m transect. We found a weak positive relationship of δ¹⁸O with HWT, which suggests that HWT cannot entirely explain species-specific differences in ¹⁸O enrichment. Instead, this can be attributed to lower water retention (i.e. higher evaporation at the same water deficit) in S. magellanicum compared to S. fuscum (Clymo, 1973; McCarter and Price, 2014). Although species differences in ¹⁸O have been reported (Aravena and Warner, 1992; Zanazzi and Mora, 2005; Bilali et al., 2013), our study suggests that the species-specific δ¹⁸O signals may not simply be a consequence of growing at different HWT but can rather reflect distinct water retention capacity in these species.

The strong influence of δ¹⁸O_precip values and, to a much lesser extent, environmental variables related to water loss, combined with a relatively small within-site variation in δ¹⁸O_tissue values, suggest that macroclimatic drivers, such as precipitation inputs, largely determine the δ¹⁸O value of peatland moss tissue. These results are promising for the use of oxygen isotopes in large-scale palaeoecological reconstructions from peat cores (Ellis and Rochefort, 2006; Chambers et al., 2012; Daley et al., 2010), although a better understanding of O isotope fractionation within tissue components and their decay relationships would improve their utility. Moreover, the simple relationships presented here can potentially be utilized to trace changes in δ¹⁸O_precip values that mirror climate variability.

5 Conclusions

Our study provides new insights into large-scale variation in Sphagnum tissue isotopic signature and suggests that isotopic composition can be used for climatic reconstructions. We show a close link between precipitation and tissue δ¹⁸O values and conclude that variation in δ¹⁸O values are mainly driven by the macroclimate, but species differences exist. In contrast, δ¹³C values were strongly related to local microtopography, while the influence of macroclimate was negligible. As suggested in earlier studies, δ¹³C values were also weakly associated with NPP. These conclusions were most strongly supported for the cosmopolitan S. magellanicum complex and species identity should be accounted for in future carbon isotope studies to avoid spurious conclusions.