Environmental Correlates of Peatland Carbon Fluxes in a Thawing Landscape: Do Transitional Thaw Stages Matter?

Peatlands in discontinuous permafrost regions occur as a mosaic of wetland types, each with variable sensitivity to climate change. Permafrost thaw further increases the spatial heterogeneity in ecosystem structure and function in peatlands. Carbon (C) fluxes are well characterized in end-member thaw stages such as fully intact or fully thawed per-mafrost but remain unconstrained for transitional stages that cover a significant area of thawing peatlands. Furthermore, changes in the environmental correlates of C fluxes, due to thaw, are not well described – a requirement for modeling future changes to C storage of permafrost peatlands. We investigated C fluxes and their correlates in end-member and a number of transitional thaw stages in a sub-arctic peatland. Across peatland-lumped CH 4 and CO 2 flux data had significant correlations with expected correlates such as water table depth, thaw depth, temperature, photosynthetically active radiation and vascular green area. Within individual thaw states, bivariate correlations as well as multiple regressions between C flux and environmental factors changed variably with increasing thaw. The variability in directions and magnitudes of correlates reflects the range of structural conditions that could be present along a thaw gradient. These structural changes correspond to changes in C flux controls, such as temperature and moisture, and their interactions. Temperature sensitivity of CH 4 increased with increasing thaw in bivariate analyses, but lack of this trend in multiple regression analyses suggested cofounding effects of substrate or water limitation on the apparent temperature sensitivity. Our results emphasize the importance of incorporating transitional stages of thaw in landscape level C budgets and highlight that end-member or adjacent thaw stages do not adequately describe the variability in structure-function relationships present along a thaw gradient.


Introduction
Northern permafrost regions contain approximately 50 % (1672 Pg) of the world's soil carbon (C) pool and peatlands store 277 Pg of this C (Schuur et al., 2008;Tarnocai et al., 2009).Thawing permafrost is projected to act as a positive feedback to climate change via the release of soil C to the atmosphere and the magnitude of this feedback remains uncertain (Schuur et al., 2013).Peatlands in the permafrost regions are currently experiencing increased rates of thaw and related changes to abiotic and biotic components (structure) and elemental cycling (function; Camill, 2005;Osterkamp, 2005).Thawing peatlands are a mosaic of different wetland types, ranging from permanently frozen (e.g., palsa) to permafrost-free and minerotrophic fens (Luoto et al., 2004).Each component of these heterogeneous landscapes has distinct C function, contributing to uncertainties in estimating landscape level C budgets.Constraining the spatial variability in peatland C fluxes and related abiotic and biotic factors, is an essential step toward estimating the positive feedback potential of thawing permafrost on climate change.
Permafrost thaw in peatlands is associated with marked changes in ecosystem structure and function.Initial ground subsidence from thaw results in wet habitats due to a high water table (Smith et al., 2012).Relative to dry areas, the seasonal frost table thaws faster in wet areas, further increasing lateral flow of water to wet areas (Quinton et al., 2009).The increase in water table depth (WTD) leads to a vegetation shift toward wetter communities and an increase in graminoid species (Camill, 1999;Camill et al., 2001;Malmer et al., 2005).Rapid changes also occur in the microbial community, notably an increased activity of methane (CH 4 ) and nitrogen cycling (Mackelprang et al., 2011).Associated with structural shifts, several functional changes Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Malhotra and N. T. Roulet: Environmental correlates of peatland carbon fluxes have been observed in thawed permafrost peatlands.Typically plant productivity increases (Vogel et al., 2009), but so does autotrophic and heterotrophic respiration (Hicks Pries et al., 2013).Organic matter decomposition may decrease due to increased anoxic conditions after ground subsidence (Camill et al., 2001).Thus, there could be an initial increase in organic matter accumulation as a result of permafrost thaw (Turetsky et al., 2007;Vitt et al., 2000).Subsequently, decomposition may increase due to the increase in easily decomposable litter from community shifts toward more vascular plants (Hodgkins et al., 2014;Turetsky, 2004) and quantity of litter (Malmer et al., 2005).Regardless of the initial increase in C accumulation, the net radiative forcing of a recently thawed area is offset by an increase in CH 4 emissions (Johansson et al., 2006;Sitch et al., 2007;Turetsky et al., 2007).This increase in CH 4 emissions may be a direct result of increased temperature on microbial processes or indirect consequences such as increases in plant mediated transport of CH 4 by increased graminoid abundance and increased anaerobic decomposition due to a high water table (Christensen et al., 2004).In addition to magnitude, the dominant pathway of CH 4 production is also altered after thaw, shifting from CO 2 reduction (hydrogenotrophic) to acetate cleavage .(acetoclastic;Hodgkins et al., 2014;McCalley et al., 2014).The change in pathway is likely related to shifts in vegetation, for example, a decrease in Sphagnum abundance could lead to an increase in pH and related increase in acetoclastic methanogens (Hines et al., 2008;Ye et al., 2012).Dissolved organic matter (DOM) is also more labile in the more thawed stages, and there is increased export of DOM out of the peatland catchment (Hodgkins et al., 2014;Olefeldt andRoulet, 2012, 2014).Recent studies highlight the interactive controls on C fluxes, emphasizing that net radiative forcing of a thawing ecosystem depends on non-linear interactions among temperature, degree of anoxia and organic matter chemistry (Lee et al., 2012;Treat et al., 2014).For example, while temperature sensitivity of CH 4 flux increases in wet habitats (Olefeldt et al., 2013), ecosystem respiration is more sensitive in dry conditions (McConnell et al., 2013).Interactive controls on C fluxes are further complicated by the variable structural conditions as thaw progresses, and the overall effect on landscape level fluxes remains unconstrained.
Changes to peatland structure and function due to permafrost degradation have been studied using a chronosequence approach with sites that have intact permafrost, completely thawed permafrost and one or two intermediate stages (e.g., Bäckstrand et al., 2010;Turetsky et al., 2007;Vogel et al., 2009).While end-member and major thaw stages of the permafrost gradient have been well characterized for plant community structure and carbon cycling, the same is not true for the transitional vegetation communities.Carbon cycling in thawing permafrost regions is spatially heterogeneous (e.g., Belshe et al., 2012;Morrissey and Livingston, 1992;Zhang et al., 2012) and a significant portion of the landscape is in varying stages of thaw.Spatial heterogeneity and transitional stages are therefore important to the ecosystem level C exchanges.It is unclear whether 3 to 5 thaw classes of intact, intermediate and fully thawed permafrost can be used to adequately extrapolate landscape scale C fluxes and their abiotic and biotic correlates.Additional thaw stages may help resolve landscape scale C fluxes in models.
Our study aims to identify the abiotic and biotic factors (hereafter, correlates) that relate to C function and investigate how these correlates change along end-member and transitional permafrost thaw stages.Our research questions are the following: (1) which correlates best explain the CH 4 and CO 2 fluxes across all thaw stages at a peatland where permafrost is thawing?(2) How does the importance of these correlates change along a gradient of increasing thaw?Our selection of measured correlates was based on current understanding of C flux relationships with temperature, moisture, pH, nutrients and plant biomass.Given the interactive nature of controls on C fluxes and variable structural changes after thaw, we expected to see no relationship between dominant correlates of C flux and degree of thaw.

Study site
The study site, Stordalen mire is located 10 km east of Abisko in Sweden (68 • 22 N, 19 • 03 E).The Stordalen peatland complex consists of several landscape units and wetland types.Permafrost is present in the dry hummocky sections of the peatland (palsa mire).Also present are areas where permafrost is thawing or has disappeared with vegetation communities that have been classified as semi-wet, wet, and tall graminoid (Johansson et al., 2006;Kvillner and Sonesson, 1980).Generally, the drier areas of the peatland complex are composed of species such as Empetrum hermaphroditum, Betula nana, Rubus chamaemorus, Eriophorum vaginatum, Dicranum elongatum and Sphagnum fuscum.The wetter areas consist of species such as E. vaginatum, Carex rotundata, S. balticum, E. angustifolium, C. rostrata, S. lindbergii, andS. riparium. The long term (1912-2003) mean precipitation measured at the Abisko Scientific Research Station (10 km from the site) is 303.3 mm, of which 150 mm occurs between June and September.The long term  mean annual temperature at the site is −0.5 • C but has surpassed 0 • C in the recent decades (summarized in Olefeldt and Roulet 2012, from observations made at Abisko Scientific Research Station).Smoothed mean annual temperature trends suggest a 2.5 • C increase between 1913 and 2006.(Callaghan et al., 2010).

Vegetation community and thaw stage selection
A total of 10 vegetation communities were selected across Stordalen to represent major stages along the thaw gradient.Selection of communities was based on an across site survey of dominant vegetation communities, coupled with characterization of water table depth and active layer thickness.The sequence of the 10 thaw stages was based on a survey of spring thaw depth and previously established vegetation community relationships with permafrost thaw (Johansson et al., 2006;Kvillner and Sonesson, 1980).

Gas flux measurements
Within each of the 10 selected communities, 3 collars of 0.05 m 2 area each were inserted in the peat surface and served as a seal for the manual gas flux measurements, with the exception of 2 communities that had only 2 collars each as they represent a small area in the mire (Table 1).Each community also had a PVC dip well installed to measure the water table depth.Methane flux was measured using opaque chambers of volume 9 or 18 L. Five headspace gas samples of 20 mL each were collected every 5 min over 20 min.Prior to collecting each sample, the headspace was mixed using a syringe.The collected headspace samples were analyzed within 24 h for concentrations of CH 4 using a Shimadzu GC-2014 gas chromatograph with a flame ionization detector, after separation on a HayeSep-Q packed column at the Abisko Scientific Research Station.Helium was used as a carrier gas at the flow rate of 30 mL min −1 .Injector, column and detector temperatures were 120, 40 and 120 • C, respectively.A 10-repetition run of known CH 4 standard (2 ppm concentration) was used to calibrate the GC before and after each sample run.Accu-racy of the analysis (calculated with the standard deviation of the 10 standard replicates) was ±0.1 to 0.75 %.Flux rates were then calculated using the slope of the linear relationship between gas concentrations and time.Only the relationships with a significant (p < 0.05) R 2 above 0.85 for the five time points were kept to calculate fluxes.If one of the five samples deviated from the linear fit, flux was calculated without it as long as the R 2 was greater than 0.95.Methane was measured on 7 days and 12 days in the 2012 and 2013 growing seasons, respectively.
For carbon dioxide flux measurements on the 28 collars, we used clear cylindrical polycarbonate chambers (13 L volume).The air enclosed within the chambers was mixed by fans and circulated through an infrared gas analyzer (PP Systems, Model EGM-4) that measured changes in CO 2 over 3-min measurement intervals (recording every 10 s for the first minute, and then every 30 s for the last 2 min).Over the 3 min measurement period, on average, temperature in the chamber only increased by 1.9 • C. Measurements were performed for full sun, with a mesh cover and finally with a black shroud, so that data from varying light intensities could be collected.Photosynthetically active radiation (PAR) was measured (Model LI-190SA, LI-COR ® , NE, USA) within each chamber over the sample interval.Fluxes were calculated using a linear regression of CO 2 concentration change over time.CO 2 was sampled on 10 days during the 2013 growing season.

Ancillary measurements
Each flux measurement of CO 2 or CH 4 was coupled with simultaneous measurements of soil temperature at 10 cm, air temperature, thaw depth and water A. Malhotra and N. T. Roulet: Environmental correlates of peatland carbon fluxes level) of each collar was measured using a real-time kinematic (RTK)-GPS.Vegetation composition for vascular plants was surveyed once every growing season in each of the collars recording the percent cover of each species.In 2013, vascular green area (VGA) was also measured on 4 days during the growing season using species specific formulae based on leafgeometry (Lai, 2012;Wilson et al., 2006).For each collar, the total number of green leaves per species was recorded along with width and length of 10 leaves for each species.The seasonality of VGA was modeled using a Gaussian fit and combined with a quadratic fit with elevation to extrapolate a spatially and temporally higher resolution data set for VGA.Throughout the manuscript we only use the modeled VGA.
Surface water was sampled near the collars on each CH 4 sampling day in thaw stages that had persistent water table throughout the growing season (Thaw stages 5, 7, 8, 9 and 10).Surface water samples were analyzed for pH and conductivity (Oakton ® portable pH conductivity meter) and reduced conductivity was calculated by removing H ion concentrations from the conductivity.Subsequently, samples were filtered using Whatman ® glass fiber filters (0.45 µm pore size) and analyzed for dissolved organic carbon and total nitrogen using a Shimadzu TOC-V Series Analyzer.

Data analyses
CH 4 flux -each flux measurement was log 10 transformed after adding 12 mg CH 4 m −2 d −1 to the original value (to account for all the negative fluxes).Log 10 transformation decreased the skew in the raw data and improved the linear relationship between CH 4 and other variables, allowing for the use of multiple linear regressions.Bivariate relationships with abiotic and biotic factors were explored with Spearman's rank-order correlations.To explore the relationship between environmental correlates and CH 4 flux, we used stepwise multiple linear regression.We used both additive and interactive effects to explore a best fit model, but found that interactive effects were either insignificant or had a weak contribution to the overall model.For ease of interpretation, given that our variables are already proxies for several interacting controls on CH 4 fluxes, we only included additive effects in our final model.The best fit model met the necessary assumptions of normality and homoscedasticity of model residuals.Multicollinearity was checked using variance inflation factors (VIFs), wherein any explanatory variable with VIFs greater than 2 was removed from the model.
Arrhenius plots were utilized to study the temperature sensitivity of CH 4 flux, regressing the log of CH 4 flux with inverse of temperature in Kelvin.
CO 2 flux -we combined all CO 2 flux data using nonlinear regression of a rectangular hyperbola to describe the relationship of net ecosystem exchange (NEE) and PAR (Bubier et al., 2003): The parameters are the following; -GP MAX : the maximum gross photosynthetic CO 2 capture at maximum PAR (µmol CO 2 m −2 s −1 ); Other than PAR, we expected to see non-linear relationships between CO 2 flux and WTD, thaw depth, soil temperature and VGA, but we did not find significant relationships.Instead we found linear relationships to be significant.Since our data were non-parametric, we used Spearman's correlation coefficients to quantify the link between CO 2 flux with abiotic and biotic variables.
Thaw gradient analyses -above analyses for CH 4 and CO 2 were repeated independently for each of the 10 thaw stages.Subsequently, strength (adjusted R 2 or Spearman's ρ) and direction of relationships between correlates and function variables were organized by thaw stage to observe whether there is a significant trend in changing correlates of CH 4 and CO 2 fluxes along the thaw gradient.The sequence of thaw stages along the gradient was based on a survey of spring thaw depth, as discussed in section 2.2.Multiple regressions were also performed for each thaw stage since CH 4 and CO 2 fluxes are not typically estimated using bivariate models.While the bivariate correlations identified how the dominant correlates change across the thaw gradient, multiple regressions across the thaw gradient provide a better idea of the changing interactive effects of abiotic and biotic correlates on CH 4 or CO 2 fluxes.
Lastly, we evaluated the relationship between CH 4 and CO 2 fluxes using a simple CH 4 : CO 2 flux ratio.For a standardized measure of CO 2 flux we use the GP MAX from each thaw stage.

Across peatland C fluxes and correlates
Mean and standard error of CH 4 flux across all collars from 2 years of sampling was 91.25 ± 8.17 mg CH 4 m −2 d −1 , ranging from −1.1 ± 0.3 to 370.2 ± 52.1 mg CH 4 m −2 d −1 (Table 1).Strongest bivariate relationships between CH 4 flux and abiotic variables were with elevation, water table depth, pH, VGA, thaw depth and surface water C : N (Fig. 1).Significant but weaker relationships were also found with soil temperature and Julian day.Total carbon (TC), total nitrogen (TN), conductivity and reduced conductivity did not have a significant relationship with CH 4 flux.
The best fit multiple regression model for CH 4 fluxes across the peatland included elevation, thaw depth, VGA and soil temperature, in decreasing order of contribution to the overall model, and these variables were able to explain 73 % of the variance in CH 4 flux (Table 2).An alternative model that excluded elevation wherein the adjusted R 2 drops to 0.62, is also reported as it better isolated the effects of VGA, soil temperature and thaw depth.The contribution (beta weights reported in brackets) of soil temperature (0.16) and thaw depth (−0.27) are similar in the model with or without elevation.The contribution of VGA increases from 0.26 to 0.58 when elevation is removed from the model.
Water table depth and thaw depth showed weak relationships with NEE and gross primary production (GPP) (calculated using NEE minus R eco ; Table 3).Soil temperature was also related to R eco and NEE.NEE and R eco were most strongly related to the mean growing season VGA (Table 3).

Correlate-flux relationships within the thaw stages
Along the thaw gradient, the strength and direction of bivariate relationships among environmental variables and CH 4 flux changed variably (Fig. 3).No significant trend along the thaw gradient was observed for the relationship between CH 4 www.biogeosciences.net/12/3119/2015/Biogeosciences, 12, 3119-3130, 2015  flux and elevation, VGA (vascular green area) and WTD.Significant trends were observed with the water chemistry variables of pH, C : N, TC, TN and conductivity and strength of correlations between correlate and CH 4 flux increased as the permafrost thawed.However, these data were only available for thaw stages with a water table (thaw stages 5 to 10).Soil temperature, thaw depth and Julian day, with data available for each of the 10 thaw stages, showed significant trends along the thaw gradient in their correlation with the CH 4 flux (Fig. 4a, b, and c).There was an increase in the amount of variance explained in the CH 4 flux by temperature as well as the slope of this relationship (Fig. 5).
The parameter estimates from the best fit model of across peatland-lumped flux data, as shown in Table 2, were used as inputs for individual multiple regression models for each thaw stage.The interactive effects of elevation, soil temperature, thaw depth and vascular green area (VGA) showed varying results across the thaw gradient (Fig. 6).The model R 2 values ranged from 0.09 (insignificant) to 0.79 (significant with p < 0.0001).Generally, elevation, soil temperature, thaw depth and VGA were better predictors of variance in CH 4 fluxes in the later stages of thaw.Model fit was nonsignificant for stages 2 and 3, and therefore their slope coefficients are not reported in Fig. 6.
The relationship of NEE with temperature and with PAR varied across the thaw gradient without a statistically significant trend.Generally there is a trend of increasing, GP MAX , α and A, going from less thawed to more thawed stages (Table 4).Furthermore, the amount of variance of NEE explained by PAR was typically higher in the more thawed stages.

Relationship between NEE and CH 4
GP MAX and CH 4 were positively correlated (ρ = 0.56, p = 0.0021; Fig. 7).We found that the best explanatory variables for CH 4 : GP MAX ratio were Sphagnum percent cover (ρ = −0.72,p = 0.008; and graminoid VGA (ρ = 0.63, p = 0.0004).While graminoid VGA did not have any interactive effects with abiotic variables in explaining CH 4 : GP MAX ratio, Sphagnum cover and soil temperature had a significant interactive effect (Table 5)

Discussion
We identified the major abiotic and biotic correlates of the ecosystem -atmospheric exchanges of CO 2 and CH 4 across Stordalen mire and found that, as per our expectation, these environment-function relationships changed variably across the thaw gradient, suggesting that correlates of CO 2 and CH 4 fluxes in transitional stages are not necessarily represented well by correlates of the end-member or adjacent thaw stages.Contrary to our expectation, we did see significant trends with thaw in certain bivariate correlations of CH 4 fluxes such as temperature sensitivity, seasonality and effect of deepening frost table during the growing season.However, these trends were absent when multiple correlates were considered together, suggesting that dominant controls on C fluxes and their interactions, change variably as thaw progresses.

Across peatland correlates of C fluxes
Strongest environmental factors associated with the CH 4 flux across all sampled collars were-elevation, water table depth, pH, VGA, thaw depth and surface water C : N.Each of these correlates is a possible proxy of one or more controls of temperature, moisture and substrate quantity and quality on CH 4 flux.Our correlations support previous findings from various wetland types (as reviewed in Lai, 2009;Olefeldt et al., 2013;Turetsky et al., 2014).The multiple regression analysis of lumped data across Stordalen also showed similar trends to other temperate, boreal or arctic peatlands.For example, a peatland complex sampled in the discontinuous permafrost region of Manitoba, Canada by Bubier et al. (1995) showed a best fit model including WTD, water chemistry and vegetation variables explaining 81 % of the variance in CH 4 fluxes.Bubier et al. (1995) reported WTD as being the strongest individual correlate with the CH 4 fluxes, but in our best fit model, WTD was not an important variable likely because the stages with little or no thaw had no water table.Elevation seems to be a better proxy for soil moisture (and other CH 4 controls), showing the highest contribution to the best fit model (Table 2).Elevation has been previously recognized as an integrator of multiple structural changes resulting from permafrost thaw and is a potentially useful component of models estimating C flux in permafrost landscapes (Lee et al., 2011).Rerunning the best fit model without elevation decreases the overall model fit by 10 % but increases the contribution of VGA to the model, while the contribution of thaw depth and soil temperature remain the same.Removal of elevation from the model better isolates the relative effects of thaw depth, temperature and VGA on CH 4 fluxes and suggests that the strongest contribution is from VGA, followed by thaw depth and soil temperature.VGA is likely a strong effect as it is linked with spatial and seasonal changes in substrate availability, litter input and root exudates and thus relates to both spatial and temporal variability in CH 4 flux (Whiting and Chanton, 1993).Soil temperature and thaw depth are significant variables in our multiple regression model of CH 4 flux, while Julian day is not.It may be that across the permafrost gradient, Julian day does not capture seasonality as well as a combination of thaw depth and soil temperature (Table 2), suggesting the role of the variable seasonal trajectories of thaw depth and soil temperature in the different thaw stages, for predicting CH 4 flux.
As expected, PAR was the strongest correlate of NEE both in across peatland-lumped data and within thaw stage data.Using the rectangular hyperbola fit, Frolking et al. (1998) reported parameters from 13 peatland sites wherein PAR explained, on average, 68 % (ranging from 47 to 89 %) of the variance in NEE.Comparatively our across peatland-lumped data fit to the rectangular hyperbola model explain a lower percent of the variance (52 %) in NEE, likely due to biases introduced by the high spatial heterogeneity on our site (Laine et al., 2009).WTD, thaw depth and soil temperature also show significant but weak relationships with CO 2 fluxes.Vascular green area seems to be a better proxy than WTD, thaw depth and soil temperature for controls on CO 2 fluxes, which makes sense as VGA represents the amount of photosynthesizing area as well as approximates above and belowground biomass which is related to autotrophic and het- erotrophic respiration (Schneider et al., 2011;Wilson et al., 2006).

Trends with increasing thaw
Bivariate relationships between correlates and CH 4 flux progress variably as the permafrost thaws, although some significant trends of increasing correlations are seen in soil temperature, thaw depth and Julian day, as thaw progresses (Figs. 3, 4 and 5).We found that as the permafrost thaws, temperature sensitivity increases (Fig. 4a), increasing thaw depth has an increasing effect on CH 4 fluxes (Fig. 4b), and there is a stronger seasonality effect (Fig. 4c).This trend of increasing correlation could be partly due to the increasing magnitude and variance of not only CH 4 fluxes but also the environmental variables with thawing permafrost.Additionally, higher VGA later in the growing season could also be result in a stronger seasonality effect (Fig. 4c) in the later stages of thaw, especially as these stages had the highest sedge VGA.The Arrhenius plots of soil temperature and CH 4 fluxes showed increased temperature sensitivity from less thawed to more thawed stages, with the slope and R 2 of this regression increasing (Fig. 5).Changing temperature sensitivity in our results contradicts results from Yvon-Durocher et al. ( 2014) that suggest a consistent temperature sensitivity of CH 4 fluxes across scales.Apparent temperature sensitivity can be confounded due to changes in substrate availability (Kirschbaum, 2006).Increasing temperature sensitivity with thaw in our results could be related to higher substrate availability (supported by higher VGA) in thawed stages switching CH 4 production from being substrate limited to becoming temperature limited.Lower temperature sensitivity in the intact permafrost could also be related to DOC quality.Olefeldt et al. ( 2012) report higher aromaticity in DOC exported from palsa and bog catchments at Stordalen compared to fen catchments and a high proportion of aromatic compounds in litter is generally associated with decreased temperature sensitivity (e.g., Erhagen et al., 2013).High temperature sensitivity in wetter sites has also been reported by Olefeldt et al. (2013) in a meta-analysis of CH 4 emissions from terrestrial ecosystems worldwide.Christensen et al. (2003) found that temperature is a limiting factor only when the WTD is 10 cm or less below the surface, whereas a lower WTD is more sensitive to WTD fluctuations than to soil temperature fluctuations.This is generally supported in our results with stages 4 to 6 that have growing season mean WTD greater than 10 cm (Table 1) having lower sensitivity to WTD than stages 7 to 10, though there is variability in both classes (Fig. 3).Our estimated temperature sensitivity for each thaw stage is the net effect of temperature on methanogenesis and methanotrophy and since we only measure the net CH 4 flux we cannot isolate the relative temperature sensitivities for the two processes.Also interesting is the effect of the increasing thaw, over the growing season, on CH 4 flux-more significant in the wetter more thawed stages than the drier intact permafrost stages (Figs. 3  and 4b).A similar trend was also emphasized in Olefeldt et al. (2013).The deepening frost table is related to temperature and thus could also represent the larger temperature sensitivity of CH 4 in later thaw stages.Additionally, larger variance in thaw depths of later thaw stages could explain the Significance of the R 2 is denoted by asterisk ( * for p < 0.05, * * for p < 0.01, and * * * for p < 0.001).Soil temperature was not a statistically significant estimate for any of the thaw stages.Elevation was significant for thaw stages 1, 4, 8 and 9; thaw depth for 1, 4, 5, 7, 8 and 10; and VGA for 8 and 9.
larger effect of thaw depth in these stages .The larger variance in thaw depth could be attributed to a steeper drop in thaw depths as the growing season progresses in the wetter thaw stages due to the dependence of thermal conductivity of peat on the degree of wetness (Quinton et al., 2009).While bivariate relationships between correlates and C flux provide insight into the possible controls on these fluxes, multiple regressions better demonstrate the interactive nature of these correlates.Re-running the best fit model of the lumped data (Table 2) for each thaw stage showed that the strength of the overall model and the parameter estimates are variable along the thaw gradient (Fig. 6).While elevation had a strong effect in across peatland-lumped data, it makes sense that it was a significant effect only for a few within thaw stage analyses (thaw stages 1, 4, 8 and 9) because these stages had diverse habitats with spatially varying elevations.Soil temperature was not a statistically significant estimate for any of the thaw stages, possibly because elevation and thaw depth are better proxies for the long term thermal regime and also relate to several other controls of CH 4 flux, as previously mentioned.VGA was only a significant effect within thaw stages 8 and 9.These results emphasize that spatial differences in elevation are not as important within thaw stages as they were in across peatland-lumped data.Also, thaw depth and VGA have variable effects but generally stronger in the thawed stages.
Similar to CH 4 flux, the strength of the major correlate for CO 2 flux (PAR) changes variably across the thaw gradient.While the across peatland relationship of NEE with temperature and PAR is weaker than that found in other peatland sites (e.g., Bubier et al., 2003), when broken down into thaw stages, the percent variance of NEE explained increases (up to 91 %; Table 4) for many thaw stages.The sample size for each thaw stages is different making it problematic to statistically compare the thaw stages.However, we found that the R 2 is not significantly correlated to the sample size for that thaw stage, suggesting that there are other factors increasing the control of PAR on NEE as permafrost thaws such as increased photosynthesizing biomass (reflected by increasing VGA and GPP).If VGA is no longer limiting, PAR sensitivity could be increasing as permafrost thaws.This is supported in the parameter fitting for each thaw stage (Table 4), the general trends observed are that the CO 2 fixed at maximum PAR (GP MAX ) increases as permafrost thaws as does the amount of CO 2 fixed per unit of PAR (α), both of which could be related to increase in VGA but also the photosynthetic capacity change from plant species changes.The amount of respiration at 0 • C generally increases with thaw, which could be related to increasing substrate availability.Trends of increasing GPP and ecosystem respiration with permafrost thaw have been reported in previous studies (e.g., Dorrepaal et al., 2009;Hicks Pries et al., 2013).However, in our results these trends are not significant along the thaw gradient and progress variably.

Relationship between NEE and CH 4
NEE is thought to be related to CH 4 emissions due to the shared association with recently produced substrate availability, root exudates and turnover and litter input, and this link has been observed in several studies (Bellisario et al., 1999;Ström and Christensen, 2007;Whiting and Chanton, 1993, etc.).In our thaw stages, there was also an overall significant and positive relationship between growing season averages of GP MAX and CH 4 (ρ = 0.56, p = 0.0021; Fig. 7).Interestingly, thaw stages 8 to 10 (graminoid dominated) have a different relationship of GP MAX and CH 4 compared with thaw stages 1 to 7 (moss dominated) suggesting a shift in the partitioning of C loss from the system as CO 2 or CH 4 with increasing thaw and changing vegetation.We expected this shift to be related to an increase in graminoid VGA (increase in lability and CH 4 emission via aerenchyma), which was supported by our data.Surprisingly, we also found the shift to be related to a loss of Sphagnum cover, perhaps due to an increase in pH and decrease in organic matter lability.Furthermore, there was a significant interaction between soil temperature and Sphagnum cover in a linear model explaining CH 4 : CO 2 , suggesting that the relationship of CH 4 and CO 2 depends on Sphagnum abundance but the effect of Sphagnum varies by temperature (

Variable changes in ecosystem relationships with increasing thaw
Permafrost thaw increases magnitude and variance of CO 2 and CH 4 fluxes as well as changes the abiotic and biotic correlates of these fluxes.As a result, the relationships between the correlates and C fluxes change.While in the lumpedacross peatland data, spatially variable factors are the dominant correlates of CO 2 and CH 4 fluxes (elevation being the best proxy for thermal regime, soil moisture, VGA, etc.), within thaw stages it is the correlates with high temporal variations that play a critical role (Julian day, deepening frost table and soil temperature).The changing correlates of CO 2 and CH 4 fluxes are important to consider from a context of upscaling these processes from within thaw stage to site to landscape scales.Changing sensitivity of CH 4 fluxes to temperature, likely related to a shift from substrate to temperature limitation going from low biomass and low nutrient palsa stages to high biomass and high nutrient thawed stages.
Based on the range of temperature response curves of CH 4 flux across the thaw gradient (Fig. 4a), applying one activation energy value to estimate landscape level CH 4 fluxes at Stordalen would not be appropriate and would likely require a set of parameterizations for the various thaw stages.
Variable temperature sensitivities to C fluxes have been recognized in major thaw stages in the past (e.g., Lupascu and Wadham, 2012), but our study demonstrates that this variability is present even in the transitional stages.Furthermore, the multiple regression analyses for each thaw stages (Fig. 6) demonstrated the changing relative importance and interactive effects of dominant correlates of CH 4 flux highlighting that controls in transitional stages of permafrost thaw are not necessarily related to controls in adjacent or end-member stages.
Paleo-ecological methods were not employed to confirm the actual thaw status of the thaw stages used in our anal-yses.Rather, our space for time approach was employed to sample the major stages of thaw at Stordalen acknowledged in previous studies (Bäckstrand et al., 2010;Johansson et al., 2006;Svensson et al., 1999, etc.) -encompassing palsa (our stages 1 to 3), internal fen (our stages 4 to 6), completely thawed flow through fen (our stages 7 to 10) type habitats -while capturing the wide range of structural conditions within each of these 3 broad thaw stages.We acknowledge that structural changes due to thaw may progress variably and attempt to capture each of these pathways.For example, palsa may collapse abruptly into a wet sedge dominated habitat that then switches to a Sphagnum lawn (our thaw stage 1 progressing into stage 8 and then to stage 5).Alternatively, this progression can be gradual with a decrease in elevation of palsa (stage 1 to 2 and then 3), followed by progression into Sphagnum lawn (stage 4 and 5).Regardless, the focus of our study was the changing correlates of C fluxes along the thaw gradient and a proposed sequence of thaw stages was required to analyze these changes.

Conclusions
Our results on the environmental correlates of C fluxes interacting and changing variably with thaw suggest that using process based models or relationships between NEE and CH 4 flux to derive landscape level C fluxes would require additional information about transitional thaw stages.
Peatlands in the discontinuous permafrost zone are highly heterogeneous, especially if they are actively thawing.Our research highlights the variability observed in structurefunction relationships with permafrost thaw.Additionally, by identifying across peatland structure-function relationships that are maintained across the heterogeneous landscape our results will assist in improving regional estimates of the carbon balance and provide insight into the level of aggregation or disaggregation needed in models to capture ecosystem level response to change.

Figure 1 .
Figure1.Relationship between CH 4 fluxes and environmental variables (across peatland-lumped data).Flux data lumped from 28 sampled collars from 2 years (19 days of sampling) showed significant relationships with Julian day, elevation (above sea level) of collar, soil temperature at 10 cm depth below surface, modeled VGA (vascular green area), water table depth (WTD), pH, thaw depth and C : N of surface water.Spearman's ρ of each correlation are shown on each graph (p < 0.0001).

Figure 3 .
Figure 3. Correlation coefficients of CH 4 flux with biotic and abiotic variables within thaw stages.Changing strength and directions of correlations are shown for thaw stages 1 to 10; where stage 1 is intact permafrost and stage 10 is fully thawed.Each data point represents a correlation analysis of n = 19 days.The missing data points in WTD, conductivity, reduced conductivity, pH, TC, TN and C : N (total carbon, nitrogen and C : N in surface water) are from the thaw stages that did not have a water table or had a correlation coefficient of 0.

Figure 5 .
Figure 4. (a) Arrhenius plots for each thaw stage.(b) Linear fit between thaw depth and CH 4 flux with progressing thaw.(c) Increase in seasonality of CH 4 flux as the permafrost thaws.

Figure 6 .
Figure 6.Multiple regression of CH 4 fluxes with elevation, thaw depth, soil temperature (temp) and vascular green area (VGA) for each thaw stage (thaw stage 1 = intact permafrost, 10 = completely thawed).Model fit R 2 values are reported along with model estimates (stacked bars).Significance of the R 2 is denoted by asterisk ( * for p < 0.05, * * for p < 0.01, and * * * for p < 0.001).Soil temperature was not a statistically significant estimate for any of the thaw stages.Elevation was significant for thaw stages 1, 4, 8 and 9; thaw depth for 1, 4, 5, 7, 8 and 10; and VGA for 8 and 9.

Figure 7 .
Figure 7. Relationship between growing season means of CH 4 flux and GP MAX (gross photosynthetic CO 2 capture at maximum PAR) across all thaw stages.Thaw stages 1 (intact permafrost) to 10 (thawed permafrost) are shown using blue to red colors, respectively.

Table 1 .
Details of each thaw stage -habitat description and dominant vegetation, vascular plant green area (VGA), presence or absence of permafrost, mean growing season water table depth (from 19 days of sampling over 2 years), mean surface water pH, reduced conductivity (Cond) and C : N ratio, growing season mean soil temperature at 10 cm depth below surface and mean growing season CH 4 fluxes.Details on CO 2 flux for each thaw stage can be found in Table4.WTD was only reported for stages with a water table on more than 5 sampled days.All values reported after ± are standard errors.

Table 3 .
Abiotic and biotic relationships with gross primary production (GPP = NEE − R eco ), net ecosystem exchange (NEE) and ecosystem respiration (R eco ).R eco has negative values and therefore, negative correlation signifies that larger VGA or soil temp have higher R eco .Similarly, since thaw depth has negative values, negative correlations with GPP and NEE mean that deeper frost tables relate to greater GPP and NEE.

Table 4 .
Parameter estimates ± standard error from rectangular hyperbola fit of NEE and PAR.PAR dependence of NEE generally increases and becomes more significant from less to more thawed stages as shown by the adjusted R 2 (p < 0.0001 for all thaw stages).There is a general trend of increasing, GP MAX , α and A, going from less thawed to more thawed stages.

Table 5 .
Multiple regression model of CH 4 : CO 2 explained by percent cover of Sphagnum and soil temperature.(R 2 0.82, F 3,11 = 11.9, p = 0.002).Slope estimates and standard errors (SE) are reported along with t ratio and associated p value.