Plant belowground traits indicate increased plant-mediated methane transport along a peatland permafrost thaw gradient

Permafrost thaw alters methane (CH₄) fluxes from subarctic peatland ecosystems. Collapsing permafrost palsas change hydrology, interstitial oxygen availability, and vegetation composition, and each of these factors contribute to net CH₄ flux by influencing CH₄ production, consumption and transport. Changes in plant-mediated CH₄ fluxes have mostly been estimated using aboveground characteristics, such as biomass and leaf area, leaving belowground parts (roots and rhizomes) understudied despite their direct contact to depth-dependent CH₄ flux processes. Here, we explored the potential of using root and rhizome traits as proxies for plant-mediated CH₄ cycling along a peatland permafrost thaw gradient in subarctic Sweden. We investigated changes in root and rhizome biomass, surface area (SA), diameter, tissue density (TD), and specific root length (SRL) along the permafrost thaw gradient, and how these traits relate to early-, middle-, peak- and season median CH₄ fluxes. We utilized chamber CH₄ flux and pore water CH₄ concentration and isotopic measurements during the productive season. Shrub SRL, diameter and isotopic data suggested increased plant-mediated carbon substrates available for acetoclastic methanogenesis across the thaw gradient. Root TD, proxy for root porosity, decreased with thaw and had negative correlations with CH₄ fluxes throughout the season. Simultaneously, herbaceous rhizome SA-CH₄ flux correlations were positive and pore water CH₄ concentrations were lowest in the fully thawed stage. These results indicated increasing herbaceous plant-mediated transport of acetoclastically-produced CH₄ with thaw. Our study demonstrates that integrating plant belowground traits with environmental and biogeochemical data can help improve CH₄ flux predictions in thawing landscapes and revealed key mechanistic insights regarding the interplay between substrate availability for methanogenesis and gas transport.

1 Introduction

Northern peatlands are significant net methane (CH₄) emitters primarily due to the predominance of anoxia throughout the peat depth profile (Gorham, 1991; McGuire et al., 2009; Zhuang et al., 2004). Owing to high organic carbon accumulation, high-latitude peatlands store approximately a third of the global carbon stock (Gorham, 1991), with a substantial fraction preserved in and beneath permafrost (Schuur et al., 2008). In the high latitudes, increasing air temperatures drive permafrost thaw, which can expose the previously frozen organic matter to soil decomposers, including methanogens. Together with altered hydrological and thermal conditions, these processes may increase CH₄ emissions and lead to positive climate feedback (Olefeldt et al., 2013; Schuur et al., 2008). Palsa mires with permafrost plateaus (palsas) are especially vulnerable to increasing temperatures and changes in precipitation (Luoto et al., 2004). Palsas are approximately 0.5–10 m high peat hummocks with a frozen core of peat and mineral soil. They form as a result of water-saturated peat freezing in winter and the dry surface peat protecting the frozen peat in summer, raising the peat above the surrounding mire (Luoto et al., 2004). In northern Fennoscandia, permafrost thaw has formed transitional gradients from frozen palsas and peat plateau to permafrost-free bogs and fens (Luoto et al., 2004; Olvmo et al., 2020). Following palsa collapse, soil moisture increases, shifting the system from an ombrotrophic bog-like system to a minerotrophic fen. Vegetation composition shifts from shrub-dominated on palsas to graminoid-dominated in fens, and soil redox and pH change, influencing methanogenic CH₄ production, methanotrophic CH₄ consumption, and CH₄ transport from the soil to the atmosphere, and net CH₄ fluxes (Perryman et al., 2020; Turetsky et al., 2014; Varner et al., 2022). Almost all arctic permafrost has been estimated to degrade and result in 120 ± 85 Gt of carbon emissions by 2100 (Lawrence et al., 2008; Schaefer et al., 2014). Together with the inevitable disappearance of palsa mires (Luoto et al., 2004; Olvmo et al., 2020), this highlights the urgent need to investigate the network of CH₄ flux controls in thawing palsa mires.

One of the key biotic CH₄ flux drivers in permafrost peatlands are wetland plants and the traits reflecting their adaptive and resource acquisition strategies. In general, plants contribute to CH₄ fluxes by mediating gas transport, providing carbon substrates for methanogens, and modifying the physicochemical properties of the rhizosphere (Chanton and Dacey, 1991; Laanbroek, 2010; Määttä and Malhotra, 2024; Noyce et al., 2014). As an adaptive strategy to anoxia, many wetland plant species develop aerenchyma which are air spaces within the plant stem, roots and rhizomes that transport oxygen (O₂) from the atmosphere to the roots and the surrounding soil (Armstrong et al., 1991; Končalová, 1990; Vroom et al., 2022). This plant-mediated gas transport mechanism enhances rhizosphere oxygenation and CH₄ consumption, but also simultaneously increases CH₄ transport from anoxic peat to the atmosphere (Laanbroek, 2010; Noyce et al., 2014; Ström et al., 2005). While peatlands store vast quantities of carbon, much of it is too recalcitrant to be directly utilized by methanogens particularly in deep peat (Hogg, 1993; Ström et al., 2005). The release of labile carbon substrates such as acetate and other organic acids by plant roots stimulates CH₄ production and emission (Joabsson and Christensen, 2001; Ström et al., 2003).

Studies using plant functional traits to predict CH₄ fluxes have primarily focused on aboveground features, such as plant green area, even though plant belowground biomass can exceed that aboveground especially in (sub)arctic ecosystems (Iversen et al., 2015; Malhotra and Roulet, 2015). Often, the aboveground traits have been used to estimate belowground characteristics, such as root porosity, based on the plant species or functional type (PFT) (e.g., Bouchard et al., 2007; Greenup et al., 2000; Iversen et al., 2015). However, above- and belowground plant traits and the functional connections between them have not been adequately studied in relation to carbon cycling and the use of aboveground traits as proxies for plant-mediated belowground carbon functioning remains unresolved (Iversen et al., 2015; Määttä and Malhotra, 2024).

Many studies have used above- and/or belowground biomass as CH₄ flux predictors but the direction of the relationship varies strongly between studies, species, and PFTs (e.g., Greenup et al., 2000; Joabsson and Christensen, 2001; Koelbener et al., 2010). While above and belowground biomass may reflect plant productivity, quantification of plant-mediated CH₄ fluxes may be better served by a knowledge of belowground traits more directly related to CH₄ production, consumption and transport. For example, root exudation in some instances can explain CH₄ flux variation better than biomass (Ström et al., 2005; Waldo et al., 2019). Because belowground traits, such as rooting depth, change with a deepening active layer (Blume-Werry et al., 2016; Iversen et al., 2015), assessing these trait variations and their effects on CH₄ fluxes along permafrost thaw gradients is essential.

Wetland CH₄ fluxes can be driven by plant-mediated carbon release and CH₄ and O₂ transport (Knox et al., 2021), which can be reflected in root and rhizome traits. Increases in root and rhizome tissue porosity (i.e., lower tissue density and larger diameter; Freschet et al., 2021a; Visser et al., 2000) can enhance plant-mediated gas transport and, depending on species and environmental conditions (e.g., water table level), can lead to increased CH₄ transport and net flux (Bhullar et al., 2013; Henneberg et al., 2016; van den Berg et al., 2020) or rhizospheric CH₄ oxidation and decreased net CH₄ flux (e.g., Fritz et al., 2011; Noyce et al., 2023). CH₄ production can be stimulated by increased labile carbon provision particularly from root exudation (Ström et al., 2003, 2012; Waldo et al., 2019) but also from rhizome and root decomposition (Hackney and De La Cruz, 1980; Scheffer and Aerts, 2000). Specific root length (SRL), the relationship between root length and dry mass, could be a proxy for carbon substrate provision, as high SRL is associated with low tissue density, shorter lifespan, high turnover rates (Bergmann et al., 2020; Eissenstat et al., 2000), and root exudation (Wen et al., 2022). However, studies on the relationships between different root traits, especially SRL and tissue density, in wetlands are rare (Iversen et al., 2015; Pan et al., 2019).

Here, we investigated (1) how root (biomass, surface area, diameter, tissue density, and SRL) and rhizome (biomass, surface area, diameter, and tissue density) traits vary along a peatland permafrost thaw gradient, (2) the correlation of these belowground traits with CH₄ fluxes, and (3) whether root and rhizome traits are strong predictors of CH₄ flux during early, middle and peak CH₄ flux during the productive season. We expected that root biomass, surface area and SRL would be highest at the intact and partly thawed permafrost stages due to lower soil moisture and high species diversity (Hough et al., 2022; Iversen et al., 2015; Keuper et al., 2017). We hypothesized that shrub rhizomes would be most abundant in the intact permafrost stage, whereas most herbaceous rhizome biomass would be found in the more thawed permafrost stage (Hough et al., 2022). We hypothesized that root and rhizome tissue density would decrease from intact to fully thawed permafrost due to increased soil moisture. We also expected that SRL would increase CH₄ flux (Saarnio et al., 2004; Ström and Christensen, 2007; Wen et al., 2022), while tissue density would have negative relationships with CH₄ fluxes throughout the productive season. We expected the trait-CH₄ flux relationships to be stronger for herbaceous plants than for shrubs at the middle and peak CH₄ flux periods, because herbaceous plants increase root biomass throughout the growing season in the tundra (Billings et al., 1978; Kummerow and Russell, 1980; Wang et al., 2016) whereas shrub root biomass can remain relatively constant (Wang et al., 2016).

2 Materials and methods

2.1 Study area

The study area is the Stordalen mire in northern Sweden (68°22^′ N 19°03^′ E; Fig. 1). The mean annual temperature of the area was −0.5 °C and mean annual precipitation sum was 204 mm for period 1985–2023, based on measurements from a nearby meteorological station (Polarforskningssekretariatet and SMHI, 2025a, b). The area contains a permafrost thaw gradient ranging from palsa with intact permafrost to partly and fully thawed permafrost, each of which includes hydrological conditions and vegetation compositions typical for each thaw stage. The permafrost is largely intact in drained palsa with ericaceous shrubs (e.g., Empetrum nigrum, Andromeda polifolia, Rubus chamaemorus), hummock-forming Eriophorum vaginatum, mosses (e.g., Dicranum sp.) and lichens (Holmes et al., 2022; Hough et al., 2022; Malhotra and Roulet, 2015). In the partly thawed hydrologically perched stage, the vegetation is dominated by Sphagnum spp. and feather mosses as well as sedges, such as Eriophorum vaginatum and Carex rotundata, typical for ombrotrophic bogs (Holmes et al., 2022; Malhotra and Roulet, 2015). The fully thawed stage includes wet minerotrophic fens that have generally higher pH and their dominant vegetation includes Eriophorum angustifolium, Equisetum fluviatile and Carex rostrata (Holmes et al., 2022; Johansson et al., 2006; Malhotra and Roulet, 2015). In general, moisture increases from intact to fully thawed permafrost stages, due to improved connections to the ground water (2007–2017 mean annual water table level ± standard deviation at partly thawed stage: 9.1 ± 4.6 cm below peat surface, fully thawed stage: 2.1 ± 4.4 cm above peat surface) (Crill et al., 2023; Holmes et al., 2022; Malhotra and Roulet, 2015). The peat depth reaches down to 3 m and is underlain by silt (Johansson et al., 2006).

Figure 1 Location of the Stordalen mire and CH₄ flux and peat coring plots within the study area. Based on vegetation classification in 2014, the study area was classified into intact permafrost (tall shrub and hummock, in light blue), partly thawed permafrost (semi-wet and wet, in dark gray), fully thawed permafrost (tall graminoid, in dark blue), open water (bright blue) and other (rock and anthropogenic surfaces, in white) (Palace et al., 2018, 2022; Varner et al., 2022). Rectangles indicate automated chamber CH₄ flux measurement plots, circles indicate peat coring plots, and the cross represents the eddy covariance tower (ICOS: SE-Sto).

2.2 Peat cores

We took peat cores using a push corer (d= 12 cm, length = 45 cm) to obtain plant material for belowground trait measurements. Despite destructive sampling being restricted, we were able to obtain three cores per thaw stage (n= 9 cores in total). In addition, peat coring was not possible next to the CH₄ flux chambers and within the eddy covariance footprint, so we estimated the vegetation composition (vascular plant and moss species) in the chamber collars with vegetation surveys (point intercept method; Sect. A1, Table C1) and found representative vegetation patches for each chamber outside of the footprint. In the plot selection, emphasis was put on vascular plant species and secondarily on moss species cover. The peat cores were collected on the 21 July 2023. The tightly packed samples were put into a cooler with ice packs and stored frozen for four days. The samples were transported from Sweden to University of Zurich, Switzerland (11 h), and stored in a freezer in −20 °C until further processing.

2.3 Plant green area measurements

To compare the chamber and peat coring plots, and to provide a non-destructive aboveground plant trait comparison to the belowground traits, we estimated species-specific vascular plant green area (GA; Tables 1 and C2). We randomly assigned three (peat coring plots) to five (chamber plots) subquadrats (16.8×16.8 cm) of the vegetation survey grid (Sect. A1) for GA measurements using the function sample in R (R Core Team, 2025). We measured green leaf and stem dimensions (leaf width and length, stem length and diameter) per species within each subquadrat. We measured 10 leaves per species per subquadrat, with the exception of fully thawed stage chamber plots where we measured all leaves and stems. For E. nigrum, we estimated the total number of leaves by counting all the leaves of three stems roughly representing the majority of the stems within the subquadrat. We calculated species-specific leaf and stem area (m²) based on geometric formulae, and GA was estimated based on the average number of green leaves and stems per m² multiplied by green leaf and stem area per m² (Wilson et al., 2007). Since the leaves of Rubus chamaemorus did not fit into a simple geometric formula, we estimated their green area using photos taken above the survey grid and ImageJ.JS v0.5.8 (Ouyang et al., 2019; Schneider et al., 2012) (see mean GA per species in Table C2).

Table 1 Vegetation characteristics at CH₄ flux chambers and corresponding peat coring plots at each permafrost thaw stage (intact, partly and fully thawed). “Core” refers to the plots where peat cores were taken for root and rhizome samples. “Species” includes the species encountered in the vegetation survey (Sect. A1), while vascular plant green area (GA) measurements also included species that were not included in the plant cover measurements. “Total” in GA is the mean vascular plant GA across plots per thaw stage. See vegetation coverage at the plot level in Table C1.

2.4 Environmental variables

In order to compare the abiotic environmental characteristics of the peat coring plots with the chamber plots, we measured peat temperature (°C), pH, water table level (cm) and active layer depth (cm) at the individual plots (Table 2). Peat temperature was measured 2–3 times a day (08:00 a.m., 12:00 p.m., 05:00 p.m.) at 10 cm depth over a three-day period (23–25 July 2023) using a soil temperature probe (Testo 720, SN 03679207, Testo, Titisee-Neustadt, Germany). Active layer depth was measured with a metal rod.

Table 2 Peat characteristics at CH₄ flux chambers and peat coring plots per permafrost thaw stage (intact, partly and fully thawed). “Core” refers to the plots where peat cores were taken for root and rhizome samples. See peat characteristics at the plot level in Table C1. ± denotes one standard deviation of the mean. Peat temperature was measured at 10 cm depth. WTD = water table depth (cm from peat surface), VWC = volumetric water content. Note that WTD for chamber plots is the mean annual WTD across the 2007–2017 period per thaw stage (Crill et al., 2023), and that peat moisture was estimated as overall chamber medians of the 2024–2025 period (see Sect. 2.4).

We measured pH on 25 July 2023. For the intact and partly thawed plots, we took ca. 8.5 g of peat from 0.5 m distance from the chamber collar at ca. 5 cm depth (intact) or 7–10 cm depth (partly thawed). The pH was measured using a pH meter (826 pH mobile, SN 1826001012160, Metrohm, Herisau, Switzerland) and pH electrode (Polyplast, Ref 238380, SN13208, Hamilton, Reno, USA) by mixing the peat sample with CaCl₂ solution (10 mL 0.01 mol L⁻¹ CaCl₂ per sample). At fully thawed plots, pH was measured directly from open water at ca. 5 cm depth. The electrode was always kept in liquid and when not in measurement solution, it was kept in a protective tube with KCl.

Half-hourly peat temperature data was measured with thermocouples installed at each chamber and acquired with CR10x data logger (Campbell Scientific Inc., UT, USA) (Bäckstrand et al., 2008; McCalley et al., 2014; Varner et al., 2026). As continuous soil moisture and water table level data were not available for the productive season in 2023, we used 10-year (2007–2017) water table depth (WTD) means and standard deviations per thaw stage from Crill et al. (2023) for a rough background estimate in soil moisture conditions. To estimate relationships between peat moisture and CH₄ fluxes, we also obtained soil moisture (volumetric water content, %VWC; Biasi et al., 2026) estimates for each chamber using soil moisture sensors (TMS-4 Standard, TOMST s.r.o., Prague, Czech Republic; Wild et al., 2019) at ca. 8 cm depth at ca. 1 m distance from the chambers. The soil moisture values were calibrated to peat using the myClim R package (Man et al., 2023; Wild et al., 2019). To obtain a representative estimate of the central tendency of the dominant peat moisture conditions at each chamber, we used overall chamber medians of daily-aggregated values (non-normal and distributions varied between chambers) from the May–August 2024 and 2025 periods.

2.5 Root and rhizome trait quantification

We cut the frozen peat cores horizontally into minimum 10 cm increments, and then further into four vertical parts (ca. 25 %) using a bandsaw. We randomly chose one of the quarters as a subsample for trait measurements. We thawed the subsamples in a refrigerator over 2–3 d and separated the aboveground vascular plants (see species-specific biomass in Table C3). We picked all living roots (light in color and elastic) and rhizomes (stiff, firm bark) from the subsample using tweezers and jeweler's glasses (2.5× magnification). Due to the very high amount of shrub roots, we took a further subsample (ca. 15 g fresh weight) from which we picked all shrub roots. We cleaned the roots and rhizomes with deionized water and grouped them into herbaceous and shrub PFTs based on a combination of branching patterns, color, and thickness (e.g., shrub roots were generally thinner and had higher branching density than herbaceous roots). The forb Rubus chamaemorus was grouped into shrub PFT (e.g., Blume-Werry et al., 2019), due to root identification challenges. We did not separate coarse and fine roots but instead focused on the average root system properties. The picked roots were stored in 70 % ethanol in a refrigerator (4 °C, max. 4 months).

We put the cleaned roots and rhizomes on a transparent tray filled with deionized water and scanned them with an Epson Perfection V700 Photo scanner (Epson, Japan), and analyzed them for root length, surface area and volume using WinRhizo Pro software (v. 2003b, Regent Instruments Inc., Québec, Canada). We used 1400 dpi for shrub roots with very fine roots (< 1 mm), 600 dpi for herbaceous roots, and 100–300 dpi for rhizomes. When rhizome length and diameter were not accurately calculated by WinRhizo, we measured them manually using ImageJ.JS. Then, all roots and herbaceous rhizomes were dried in 70 °C and woody shrub rhizomes in 102 °C for 4 d and weighed at 0.00001 g precision to obtain root and rhizome dry weight. To obtain a proxy of belowground gas transport potential (i.e., porosity), we calculated root and rhizome tissue density (TD, g cm⁻³) by dividing their dry weight by fresh volume. For shrub rhizomes, the TD represents all the different rhizome tissues (e.g., bark, wood and pith). While this is not a direct measurement of the gas transport ability of a woody rhizome, low tissue density could indicate increased rhizome aeration (e.g., pith cavities) and possibly gas transport (Vroom et al., 2022). Specific root length (SRL, m g⁻¹) was calculated by dividing total root length by its dry weight.

Since peat core depths varied between plots, we standardized root total length (m m⁻²), biomass (g m⁻²) and surface area (SA, m² m⁻²) to 30 cm depth by taking into account core volume and peat bulk density. We also standardized the traits to 10 cm depth for depth-specific analyses. In addition, shrub root and rhizome biomass, length and SA were scaled from the sub-subsample to the increment level and to 30 cm standardized depth. We calculated root and rhizome volumes by summing diameter-class-specific volumes from WinRhizo (Rose, 2017). All traits were calculated per PFT but also summed (biomass, length, and SA) or biomass-weighted (TD, SRL, diameter) to get an across-PFTs trait estimate. We used the rest of the peat core sample (remaining ca. 75 %) to determine peat bulk density and gravimetric water content (GWC). We thawed the sample in a refrigerator over 2–3 d, removed the aboveground vascular plants, took a 5 ± 0.2 g subsample, and took its fresh and dry weight to obtain bulk density and GWC.

2.6 CH₄ flux measurements

CH₄ flux data were obtained from automated chambers (Varner et al., 2026). The details of the chamber measurement system in general are presented in Bäckstrand et al. (2008). In 2012, the chambers were replaced to similar design as in Bubier et al. (2003). In each thaw stage (intact, partly and fully thawed), three transparent automated gas flux chambers (basal area: 0.2 m², height: 15–75 cm) measured CH₄ and CO₂ every three hours with 5–8 min enclosure time and 5 min air flushing before and after closure (Bäckstrand et al. 2008; Holmes et al. 2022). The CH₄ concentrations were measured with a Los Gatos Research Fast Greenhouse Gas Analyzer (Los Gatos Research Inc, Los Gatos, CA, USA) and the timing control and data acquisition was done with Campbell CR10x data logger (Campbell Scientific, Logan, UT, USA), which was located inside a temperature-controlled cabin. The CH₄ flux data were fit using linear regression. The data were quality-checked using regression diagnostics where CH₄ fluxes outside of significance thresholds at 95 % and 90 % confidence levels were filtered out (Snedecor and Cochran, 1989), as well as manual quality flags that were used to filter out data points affected by known instrumental issues.

We aggregated the higher frequency CH₄ flux data to daily scale by taking the daily CH₄ flux median. Median was chosen to estimate the central tendency of CH₄ fluxes due to skewed and non-normal flux data distributions, and to avoid the disproportionate influence of short-lived high CH₄ flux events that may not be related to belowground plant-mediated CH₄ cycling processes that were the main perspective of this study. Then, we subset the daily-aggregated chamber CH₄ flux data to contain the dates of the productive season in 2023. We defined the productive season as the period between the first and last passing of daily median net ecosystem CO₂ exchange (NEE) past zero, as this can be considered as the plant active season relevant for plant-mediated CH₄ fluxes (Körner et al., 2023). For NEE, we utilized the half-hourly gap-filled NEE (NEE_VUT: estimated across 40 friction velocity threshold realizations for each year according to FLUXNET data processing protocols; quality flags = 3 filtered out) from the ICOS Abisko-Stordalen Palsa Bog (SE-Sto) eddy covariance data processed with FLUXNET standards (Lundin et al., 2023; Pastorello et al., 2020). The resulting productive season used for CH₄ flux data subsetting was 19 May 2023 to 30 August 2023 (Fig. B1). As a background proxy for potential plant carbon substrate provision and plant-mediated gas transport (Knox et al., 2021), we also utilized the half-hourly gap-filled gross primary production (GPP; GPP_NT_VUT: nighttime partitioning at the daily scale and GPP_DT_VUT: daytime partitioning with light response curves at the hourly scale, with variable friction velocity thresholds implemented as an ensemble of 40 friction velocity thresholds by the FLUXNET team; see Lundin et al., 2023 and Pastorello et al., 2020 for details) data from Lundin et al. (2023). However, it is important to note that the EC tower footprint covers only intact and partly thawed stages and does not estimate the NEE and GPP for the fully thawed stage (Łakomiec et al., 2021). To explore the possible diurnal-scale plant-mediated CH₄ transport as a supporting analysis, we also aggregated CH₄ fluxes to the hourly scale using means and medians (Fig. B1).

To explore the variation in the trait-CH₄ flux relationships during the productive season, we calculated CH₄ flux period values (early, middle and peak CH₄ flux) by fitting Gaussian generalized additive models (GAM) with identity link function for daily median CH₄ flux data per chamber with function gam from package mgcv (Wood, 2011, 2017). In the models used for determining CH₄ flux period values, CH₄ flux (response variable) was transformed to approximate normality using inverse hyperbolic sine and date was included as an explanatory variable in a thin plate regression spline. Following an approach similar to TIMESAT (Eklundh and Jönsson, 2015; Jönsson and Eklundh, 2002) and Delwiche et al. (2021), CH₄ flux period values were calculated as follows: (1) peak CH₄ flux: the highest predicted CH₄ flux value during the productive season, (2) early CH₄ flux: the point where the predicted daily CH₄ flux reached 10 % of the peak CH₄ flux for the first time, and (3) middle CH₄ flux: the point where the predicted daily CH₄ flux reached 50 % of the peak CH₄ flux for the first time. We also calculated the whole season median to represent the whole productive season CH₄ flux estimate per chamber. Since one chamber at intact permafrost thaw stage had GAM R² < 0.1, we used a moving 7 d median to estimate the early, middle and peak CH₄ season values. See details of GAM fit in Table C4.

2.7 Porewater CH₄ concentration and isotopic measurements

Data of porewater CH₄ and CO₂ concentrations, δ¹³CH₄ and δ¹³CO₂ were obtained from Wilson et al. (2026). The porewater samples were collected from the field near the automated flux chambers using stainless steel sippers on 25 and 26 July 2023. The bottom of the sippers is perforated, and a syringe is attached to the top of the sipper via a stopcock. The sipper was inserted in the peat to the desired depth (0–4, 10–14, 20–24, 30–34, 40–44, 50–54, 60–64, 70–74, and 80–84 cm below peat surface) and the syringe was used to carefully draw up 30 mL of porewater. The porewater was then injected into pre-evacuated sealed glass serum vials. Phosphoric acid (1 mL, 10 %) was added to each vial to preserve the samples and ensure that all dissolved inorganic carbon was in the form of dissolved CO₂ in equilibrium with gaseous CO₂. Porewater CO₂ and CH₄ concentrations and isotopes were analyzed simultaneously via headspace analysis, corrected for solubility using Henry's Law, on a Finnigan Mat Delta V Isotope Ratio Mass Spectrometer coupled to a gas chromatograph (Merritt et al., 1995). The analytical uncertainty based on repeated measurements of a standard was < 0.15 ‰ for δ¹³CO₂ and δ¹³CH₄ and <5 % for CO₂ and CH₄ concentrations. To obtain a robust estimate of the isotopic fractionation between CO₂ and CH₄, we calculated the isotope fractionation factor (α_C) (e.g., Hodgkins et al., 2014; McCalley et al., 2014) with the following formula:

$$\begin{array}{}\text{(1)}& {\mathit{\alpha }}_{\mathrm{C}}={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{CO}}_{\mathrm{2}}+\mathrm{1000}}{{\mathit{\delta }}^{\mathrm{13}}{\mathrm{CH}}_{\mathrm{4}}+\mathrm{1000}}}\end{array}$$

2.8 Statistical analysis

To confirm the similarity of the peat core (n= 9) and CH₄ flux plots (n= 9), we used principal component analyses (function PCA from package FactoMineR; Lê et al., 2008), PERMANOVA (function adonis2 from package vegan; Oksanen et al., 2025), and Euclidean distances between plots within thaw stages, using normalized plot herbaceous and shrub plant cover and GA, mean peat temperature and pH (see Fig. B2).

For assessing general differences in daily median CH₄ flux between thaw stages across the productive season, we built linear mixed effects models where CH₄ flux was the response variable and thaw stage (factor) was the explanatory variable, and chamber (n= 3 per thaw stage) was the random effect. To account for the normality of residuals and chamber-dependent temporal autocorrelation, we applied inverse hyperbolic sine transformation to CH₄ fluxes, and modeled temporal autocorrelation with autoregressive autocorrelation structure of order 1 (corAR1). We built three different models, where each had either intact, partly or fully thawed permafrost stage as the reference level and fit the model with restricted maximum likelihood and obtained p-values (α= 0.05) from Wald t-test. Linear mixed effects models were built with function lme from package nlme (Pinheiro and Bates, 2000; Pinheiro et al., 2023).

Due to the small sample size (n= 3 per thaw stage), we were unable to test for statistical differences in belowground traits between thaw stages and sampling depths (traits, porewater CH₄ concentration and isotopic data) and relied on descriptive methods. However, to explore the general direction and strength of the relationships between root and rhizome traits and CH₄ fluxes (early, middle, peak, and season median) across the thaw stages, we used linear regressions with centered and scaled trait and CH₄ flux values. Only models which did not violate assumptions of homoscedasticity, linearity and normality of residuals were included in the analyses. In cases with little trait variation within each thaw stage, pseudoreplication was considered too severe for reliable model estimates, and these model results were not reported. We focused on model R², slopes, and model p-values. In cases where linear regressions were unreliable or violated the model assumptions, we explored the relationship strength with Kendall τ correlation coefficients. We acknowledge that the model and correlation estimates can be unstable with small sample size and thus use them as descriptive exploratory tools for estimating the changes in the general directions of the trait-CH₄ flux relationships across thaw stages during the productive season. All data processing and statistical analyses were conducted in R v4.5.1 (R Core Team, 2025).

3 Results

3.1 More herbaceous roots with decreasing tissue density

Belowground traits varied between thaw stages, and the biomass-weighted traits followed shifts from shrub- to herbaceous-dominated plant communities (Fig. 2, Table C1). Summed root biomass and surface area (SA) were 22.2× and 17.0× lower, respectively, at the fully thawed than intact stage, but herbaceous root biomass (mean: 238.3 ± 172 g m⁻²) and SA (mean: 19.9 ± 13 m² m⁻²) were highest at partly thawed stage. Total shrub root length reached a very high mean of 178.62 km m⁻² and a maximum of 199.7 km m⁻² in the intact stage, substantially higher values than for herbaceous roots (mean: 19.79 km m⁻² and max 33.33 km m⁻² in the partly thawed stage) (Fig. B4). While shrub rhizome biomass and SA decreased from intact to partly thawed stage, herbaceous rhizomes showed an opposite trend (e.g., mean SA 0.13 ± 0 to 3.41 ± 0.08 m² m⁻², 26-fold increase). Notably, herbaceous root and rhizome tissue density (TD) were 1.5× and 2.4× lower in fully thawed than in intact permafrost, respectively (mean ± standard deviation root 0.004 ± 0.0008 to 0.003 ± 0.0002 g cm⁻³; rhizome 0.2 g cm⁻³ ± 0 to 0.08 ± 0.01 g cm⁻³), and these decreases were reflected in the biomass-weighted values (root: 1.5× lower, rhizome: 7.4× lower). While specific root length (SRL) was highest at partly thawed stage (mean shrub 193.6 ± 32 m g⁻¹, 1.4× higher than at intact; mean herbaceous 71.6 ± 35 m g⁻¹, 2.3× and 1.2× higher than at intact and fully thawed, respectively), the biomass-weighted means showed a 2.2-fold decrease from intact to fully thawed stage.

Figure 2 Belowground traits varied between thaw stages. Root traits (a–e) are shown in the top row and rhizome traits (f–i) in the bottom row. Large points represent group means, small points plot-scale measurements, and error bars ±1 standard deviation of the mean. The colors represent plant functional types (PFT; dark grey: herbaceous, light grey: shrub, white: PFTs combined). Trait values for combined PFTs were summed herbaceous and shrub values (a, b, f, g) or biomass-weighted means between herbaceous and shrub PFTs (c–e, h, i).

Most of the root biomass was found in the top 10 cm for shrub (coefficient of variation CV and mean for intact: 69 %, 766.0 ± 155 g m⁻²; partly thawed: 67 %, 65.1 ± 28 g m⁻²) and herbaceous PFTs (CV and mean intact: 56 %, 17.5 ± 9 g m⁻², partly thawed: 45 %, 86.4 ± 89 g m⁻², fully thawed: 73 %, 27.7 ± 20 g m⁻²) (Fig. B3). Shrub roots reached the maximum total length of 167.3 km m⁻² at top 10 cm at the intact stage (Fig. B4). SRL and root TD varied less with depth across thaw stages (max CV SRL: 51 %, TD: 31 %; partly thawed herbaceous). All herbaceous rhizomes were present in the top 10 cm (mean biomass intact: 11.5 ± 0 g m⁻², partly thawed: 101.9±109 g m⁻², fully thawed: 107.8 ± 12 g m⁻², Fig. B5). Most of the shrub rhizomes were located in the top 10 cm (CV and mean biomass intact: 78.5 %, 392.4 ± 227 g m⁻²; partly thawed: 66.7 %, 93.9 ± 15 g m⁻²), but they were also found in 10–20 cm depth (mean intact: 109.4 ± 154 g m⁻², partly thawed: 28.9±0 g m⁻²).

3.2 Root tissue density, specific root length and diameter correlated best with CH₄ fluxes

The relationships between root traits and CH₄ fluxes (CH₄ flux difference between intact, partly and fully thawed p<0.01, fully and partly thawed p>0.05, Sect. A2) differed between PFTs but with little temporal variation. Increasing root TD was associated with decreasing CH₄ flux for herbaceous (min standardized regression beta coefficient, β at peak CH₄, max β at middle CH₄) and shrub PFTs (min β at middle CH₄, max β at CH₄ season median) (Fig. 3). Shrub root diameter had negative (min β at early CH₄, max β at CH₄ season median) and SRL positive (min β at peak CH₄, max β at CH₄ season median) associations with CH₄ fluxes. In general, the direction and strength of the relationships between root traits and CH₄ fluxes did not vary strongly across CH₄ flux season values (Figs. 3, B6).

Figure 3 Root tissue density, diameter, and specific length showed strongest relationships with CH₄ fluxes. (a–b) Relationships between root diameter and CH₄ fluxes (early, middle, peak and season) per plant functional type (PFT; a) and across PFTs (biomass-weighted means; b). (c–d) Relationships between root tissue density and CH₄ fluxes per PFT (c) and across PFTs (biomass-weighted means; d). (e–f) Relationships between specific root length and CH₄ fluxes per PFT (e) and across PFTs (biomass-weighted means; f). In (a), (c) and (e), dark gray: herbaceous plants, light gray: shrubs. In (b), (d) and (f), light blue rectangle: intact, grey circle: partly thawed, dark purple triangle: fully thawed stage. In all plots (a–f), lines represent linear regression for visualization, while R² and standardized effect size/slope (β) are based on linear regressions with centered and scaled trait and CH₄ flux values (lines, R² and β are only shown when R²>0.2, model p<0.1, and the linear regression model assumptions were met). Kendall's τ is shown for relationships where linear regression model assumptions were not met and τ≥0.2 (note that due to the low sample size, p-values should be interpreted with caution). See relationships between root biomass, surface area and CH₄ fluxes in Fig. B6.

PFTs combined, all root traits except diameter showed negative trends, where the strongest relationships were for root TD (min β at early and peak CH₄, max β at middle CH₄) (Fig. 3). In addition, summed root biomass and SA showed relatively strong negative relationships, with biomass showing slightly stronger negative relationships than SA (biomass: min Kendall's τ at early CH₄, max τ at middle CH₄; SA: min τ at early CH₄, max τ at peak CH₄) (Fig. B6).

Herbaceous and shrub rhizome traits showed contrasting relationships with CH₄ fluxes (Figs. 4, B7). Herbaceous rhizome biomass and SA showed positive relationships with CH₄ fluxes, with stronger trends for SA (SA: min β at peak CH₄, max β at CH₄ season median). Increases in herbaceous rhizome diameter were also associated with slightly higher CH₄ fluxes (min τ at early to peak CH₄, max τ at CH₄ season median). In contrast, shrubs had negative trends between rhizome biomass and CH₄ fluxes (min τ at early, middle and season median CH₄, max τ at peak CH₄; Fig. B7). While summed rhizome SA did not show clear trends (Fig. 4), higher summed rhizome biomass was associated with lower CH₄ fluxes (nonsignificant β at early CH₄, max β at peak CH₄; Fig. B7), while the relationship between rhizome diameter and CH₄ flux was positive (min τ at early and middle CH₄, max τ at peak and season median CH₄) (Fig. 4).

Figure 4 Relationships between rhizome traits and CH₄ fluxes were strongest for surface area and diameter. Columns (a–b) Relationships between rhizome surface area and CH₄ fluxes (early, middle, peak and season; rows) per plant functional type (PFT; a) and across PFTs (summed; b). Columns (c–d) Relationships between rhizome diameter and CH₄ fluxes per PFT (c) and across PFTs (biomass-weighted means; (d). In (a) and (c), dark gray color represents herbaceous plants while light gray represents shrubs. In (b) and (d), light blue rectangles represent the intact permafrost, dark blue circles represent the partly thawed stage, and dark purple triangles represent the fully thawed stage. Lines represent linear regressions for visualization, while R² and standardized effect size/slope (β) are based on linear regressions with centered and scaled trait and CH₄ flux values (lines, R² and β are only shown when R²>0.2, model p<0.1, and the linear regression model assumptions were met). Kendall's τ is shown for relationships where linear regression model assumptions were not met and τ>0.2. See relationships between rhizome biomass, tissue density and CH₄ fluxes in Fig. B7.

In comparison with belowground trait relationships, CH₄ fluxes were weakly and positively correlated with plant green area (GA) and peat temperature (Fig. B8). GA, which also did not show clear trends with permafrost thaw (Fig. B9), showed no clear relationships with CH₄ fluxes (max τ for shrub at peak CH₄; Fig. B8). Peat temperature did not show clear relationships with CH₄ fluxes (max τ at middle CH₄; Fig. B8). However, increases in peat moisture were generally associated with higher CH₄ fluxes (τ= 0.73 at early, middle and season median CH₄; Fig. B8), with the highest CH₄ fluxes occurring in the wettest plots in the fully thawed stage (median VWC = 100 %), and lowest CH₄ fluxes in the drier intact stage (median VWC = 43.4 %).

4 Discussion

Permafrost thaw in subarctic peatlands can lead to substantial alterations in peatland hydrology, vegetation composition and carbon cycling, and higher CH₄ emissions with higher water table level and more anoxic conditions belowground (Johnston et al., 2014; Knoblauch et al., 2018; Turetsky et al., 2014; Varner et al., 2022). These trends were visible in the Stordalen Mire where CH₄ fluxes increased significantly from intact to fully thawed permafrost stage with maximum CH₄ fluxes in August (Sect. A2 and Fig. B10), and porewater CH₄ concentrations increased with depth and CH₄ production slightly shifted from hydrogenotrophic to acetoclastic methanogenesis (Fig. B11), as has been observed previously at Stordalen (Hodgkins et al., 2014; Holmes et al., 2022; McCalley et al., 2014; Perryman et al., 2020; Varner et al., 2022). Shifts from hydrogenotrophic to acetoclastic methanogenesis can indicate increased root exudate-driven CH₄ production (see Sect. 4.2) (Saarnio et al., 2004; Ström et al., 2003).

The increased CH₄ fluxes with permafrost thaw in Stordalen are driven by a complex network consisting of drivers such as higher water table level and increased anoxia (Holmes et al., 2022; Malhotra and Roulet, 2015), higher peat temperatures at the end of growing season (Mollenkopf et al., 2026), increased labile carbon availability (Hodgkins et al., 2014; Holmes et al., 2022; Hough et al., 2022), shifts in microbial communities from hydrogenotrophic to acetoclastic methanogens (McCalley et al., 2014), and increasing dominance and productivity of graminoid vegetation (Christensen et al., 2004; Johansson et al., 2006; Malhotra and Roulet, 2015; Öquist and Svensson, 2002; Varner et al., 2022). For the first time in Stordalen and, to our knowledge, in other thawing permafrost peatlands, we explored the effect of belowground plant traits on CH₄ fluxes. In general, root and rhizome traits reflected plant community transitions and were consistent with an increase in plant-mediated transport of acetoclastically-produced CH₄.

4.1 Root and rhizome traits reflect plant community shifts along permafrost thaw gradient

The vegetation composition in Stordalen follows the changing hydrology and soil nutrient availability. The vegetation shifts from shrub-dominated (intact stage; e.g., R. chamaemorus, B. nana, and E. nigrum) palsas to increasingly more graminoid- and Sphagnum-dominated semi-wet bogs (partly thawed; e.g., V. oxycoccos, C. rotundata, and S. capillifolium) and wet fens (fully thawed; E. angustifolium, E. fluviatile, and S. riparium; Tables C1 and C2) (Hough et al., 2022; Johansson et al., 2006; Malhotra and Roulet, 2015; Varner et al., 2022). Following these vegetation transitions, plant traits reflect plant adaptation to the changing environmental conditions (Dong et al., 2020; Moor et al., 2017; Pan et al., 2019). While studying plant trait variation and its influence on ecosystem functioning often requires species-specificity to capture changes in diversity (Kichenin et al., 2013; Violle et al., 2014) and separation into fine and coarse root fractions (Freschet et al., 2021b; McCormack et al., 2015), CH₄ flux models often utilize PFTs which have been found to adequately represent the collective trait effects on CH₄ fluxes in peatlands (Gray et al., 2013; Laine et al., 2022). Thus, we here focused on root and rhizome traits at the PFT (herbaceous and shrub) level. This approach provides empirical constraints on plant-mediated carbon substrate provision and CH₄ transport that can inform parameterization and sensitivity analyses in process-based wetland CH₄ models, such as LPJ-GUESS (Smith et al., 2001, 2014), and HIMMELI (Raivonen et al., 2017).

Following vegetation community shifts, most belowground shrub biomass and total length were found in the intact stage while herbaceous biomass increased from intact to fully thawed stage. The shrub root biomass values were lower than previous estimates for R. chamaemorus, E. nigrum and A. polifolia (total > 2000 g m⁻²) in Stordalen (Wallén, 1986), but are much higher than recent belowground biomass estimates (∼ 280 g m⁻²; Hough et al., 2022). In addition, Hough et al. (2022) reported higher total root biomass at the fully thawed (∼ 218 g m⁻²) than partly thawed stage (∼ 146 g m⁻²), whereas this study showed the opposite (Fig. 2). These differences could arise from different root biomass standardization and sampling methods, a problem highlighted in root ecology (Freschet et al., 2021a): Wallén (1986) calculated fine root biomass based on ¹⁴C activity which can lead to overestimations (Wallén, 1986), while Hough et al. (2022) sampled a larger area (625 cm²) than in our study (∼ 79 cm²), leading to differences particularly at the fully thawed stage, where graminoids root more vertically than laterally. In addition, our root samples may have consisted of a different mix of species than in Hough et al. (2022), which together with the low sample size, may have influenced these results. Nonetheless, root biomass decreased by depth, as is common particularly for shrubs (Iversen et al., 2015; Wallén, 1986). As we observed R. chamaemorus roots up to 30 cm depth, similar to Wallén (1986) and Keuper et al. (2017) in Stordalen and by Metsävainio (1931) in Finnish peatlands, and as E. angustifolium and E. fluviatile roots have been found at 40 cm depth (Iversen et al., 2015; Metsävainio, 1931; Shaver and Billings, 1975), the maximum rooting may occur deeper than our sampling depth (ca. 30 cm). In general, these findings, and particularly the remarkably high total shrub root length values, showed that permafrost-affected peat soils are rich in plant roots and likely contribute to the carbon cycling of these ecosystems (Bardgett et al., 2014; Iversen et al., 2015).

Variation in tissue density (TD) and specific root length (SRL) (Fig. 2) with thaw suggest changes in plant gas transport and resource acquisition. Increasing soil moisture with thaw (via higher water table level; Table 2 and Fig. B8) and subsequent aerenchyma formation (e.g., Jackson and Armstrong, 1999) may have decreased TD. The increase in herbaceous rhizome SA may reflect greater aerenchyma proportion, confirming rhizome SA as a proxy for plant gas transport capacity. SRL was highest and root diameter lowest at the partly thawed stage, which could reflect increased resource acquisition particularly for shrubs (Bergmann et al., 2020), but further research with larger sample sizes is needed to investigate the potential changes in shrub resource acquisition . As the peat C $/$ N ratio is highest in the partly thawed stage (Hodgkins et al., 2014), the high SRL may reflect low N availability, shifts in species and species-specific resource acquisition (Hewitt et al., 2019; Keuper et al., 2017). However, as fine root SRL is more meaningful for resource acquisition (Bergmann et al., 2020; McCormack et al., 2015), our results based on the whole root system provide only a rough estimate. Nonetheless, 99% of total shrub root length was < 2 mm in diameter (vs. 31 %, 75 %, and 63 % in intact to fully thawed stages, respectively, for herbaceous), and thus within the traditional fine root definition (McCormack et al., 2015).

4.2 Belowground traits reflect increasing plant-mediated CH₄ transport

The increasing peat moisture and associated changes in belowground plant traits contributed to the increasing CH₄ fluxes along the permafrost thaw gradient. Given our exploratory results (see Sect. 3.2 and Fig. B8) and earlier studies conducted at Stordalen (Holmes et al., 2022; Malhotra and Roulet, 2015), the increasing peat moisture was likely one of the primary drivers of the increasing CH₄ fluxes along the thaw gradient. However, due to the low sample size we were unable to separate the effects of plant traits and peat moisture in multivariate models, and larger sample sizes are recommended to analyze the interactions between peat moisture and plant belowground traits in future studies. Nonetheless, the univariate relationships between plant belowground traits and CH₄ fluxes indicated (1) maintained carbon substrate provision and (2) increasing plant-mediated CH₄ transport across the permafrost thaw gradient. In Stordalen, permafrost thaw has led to increasing plant-derived labile carbon availability and acetoclastic methanogenesis (Hodgkins et al., 2014; Holmes et al., 2022; McCalley et al., 2014; Ström and Christensen, 2007), but plant-mediated CH₄ transport has so far been mostly assumed based on literature and plant aboveground characteristics and activity, such as number of shoots and gross photosynthetic rates (Öquist and Svensson, 2002).

Shrub roots suggested increasing carbon substrate provision. The higher CH₄ fluxes with increasing shrub SRL (and smaller diameter; Bergmann et al., 2020) may be related to increased root exudation (Guyonnet et al., 2018; Wen et al., 2022). The negative root TD-CH₄ flux relationships support this, as higher root TD may be associated with lower root exudation (Wen et al., 2022). In contrast, herbaceous SRL was weakly associated with CH₄ fluxes, possibly due to decoupling of root exudation and SRL in some Cyperaceae species in nutrient-limited environments (Wen et al., 2022), and the SRL also including coarse roots (Picon-Cochard et al., 2012). However, based on the ¹³C-depleted δ¹³CH₄ and high α_C, hydrogenotrophic methanogenesis may have dominated in the topmost peat, and most of the dissolved, less bioavailable organic carbon possibly originated from Sphagnum in the partly thawed stage (Hodgkins et al., 2014; McCalley et al., 2014; Mondav et al., 2014; Wilson et al., 2022). Nonetheless, the lower α_C, more enriched δ¹³CH₄, and higher CH₄ concentrations below the root sampling depth (30–44 cm) indicated acetoclastic methanogenesis (Fig. B11) (e.g., by Methanosarcina; McCalley et al., 2014), possibly fueled by root exudates (Saarnio et al., 2004; Ström et al., 2003; Ström and Christensen, 2007). Indeed, based on photosynthetic rates, the partly thawed stage has been estimated to have more labile carbon for methanogenesis (Öquist and Svensson, 2002), and E. angustifolium and E. vaginatum have high potential root exudation rates in Stordalen (Ström and Christensen, 2007) and graminoid root exudation rates can exceed that of shrubs in the partly and fully thawed stages (Mollenkopf et al., 2026). Based on these previous studies, it is possible that the SRL based on the whole heterogeneous root system did not represent the potential root exudation accurately for herbaceous plants. Thus, roots (including herbaceous) may have exuded labile carbon, but with shrubs having a lesser effect on the net CH₄ flux.

Herbaceous plant coverage increased with thaw, which was reflected in traits indicating plant-mediated CH₄ transport (Fig. 2, Table C1). Higher herbaceous rhizome SA was positively associated with greater CH₄ fluxes, which, together with the positive diameter-CH₄ flux trends, could indicate enhanced CH₄ transport via larger and more porous rhizomes (Määttä and Malhotra, 2024). Rhizome-mediated CH₄ transport could be particularly relevant for plants with pressurized gas transport, such as Equisetum spp., where pressure gradients between the atmosphere, leaves, and rhizomes can drive CH₄ emission via old leaves and broken stems (Vroom et al., 2022). As shrub rhizomes were not associated with CH₄ fluxes (Figs. 4, B8), they may not drive plant–mediated CH₄ cycling in Stordalen, though standardized measurements of, e.g., rhizome length and biomass are needed to confirm this (e.g., Freschet et al., 2021a).

Plant-mediated CH₄ transport may have been pronounced in the fully thawed stage with high peat moisture, herbaceous plant coverage, and low porewater CH₄ concentrations (Table C1 and Figs. B8, B11), indicating plant CH₄ uptake. As low root and rhizome TD indicates higher porosity (Ye and Ryser, 2022), the negative relationships between root TD and CH₄ flux suggest increased plant-mediated CH₄ transport in response to increasing peat moisture. The isotopic evidence (heavier δ¹³CH₄ and lower α_C) may further suggest plant-mediated transport of acetoclastically-produced CH₄ at the fully thawed stage (Chanton, 2005). This supports previous findings of higher herbaceous shoot numbers increasing CH₄ flux in the fully thawed stage (Öquist and Svensson, 2002). However, a recent study conducted in the same productive season found that 30 % of the CH₄ flux in the fully thawed stage consisted of graminoid plant-mediated CH₄ transport and the rest of CH₄ production stimulated by graminoid root exudation, while in the partly thawed stage plant-mediated CH₄ transport dominated the CH₄ flux (80 %) (Mollenkopf et al., 2026). Nonetheless, both this and the study by Mollenkopf et al. (2026) together suggest increased herbaceous plant-mediated CH₄ transport and carbon substrate provision in the fully thawed stage. In Stordalen, most of the plant-mediated CH₄ transport likely occurs via diffusion, typical for Carex spp., E. angustifolium, and E. vaginatum (e.g., Ge et al., 2025; Noyce et al., 2014; Bhullar et al., 2013). CH₄ diffusion may be more efficient through E. angustifolium with low root branching in the fully thawed than at partly thawed stage with more densely rooted Carex species (Schimel, 1995; Shaver and Billings, 1975; Table C1). However, pressurized CH₄ transport via Equisetum fluviatile (in fully thawed stage) could be negligible and requires further investigation as we did not separate the belowground material into species (Hyvönen et al., 1998; Kankaala and Bergström, 2004). As graminoid (e.g., E. angustifolium) roots have low bioavailability in the fully thawed stage (Hough et al., 2022), and root porosity and decomposability may be decoupled (Pan et al., 2019; Ye and Ryser, 2022), the low root and rhizome TD were unlikely to lead to increased decomposition and carbon substrate provision.

Simultaneously to plant-mediated CH₄ transport, herbaceous roots may have enhanced rhizospheric CH₄ oxidation. Previous studies at Stordalen have reported shifts in methanotroph communities with thaw as well as an increase in CH₄ consumption, which may result from greater plant-mediated rhizospheric oxidation (i.e., more O₂ transport with lower root TD) and acetoclastic methanogenesis (Perryman et al., 2020; Singleton et al., 2018). In addition, Sphagnum spp. can host methanotrophic CH₄-consuming bacteria (e.g., Larmola et al., 2010; Putkinen et al., 2014), which may have increased CH₄ consumption in the more O₂-rich intact and partly thawed stages dominated by S. balticum and S. capillifolium, with intact stage chambers possibly having higher CH₄ consumption rates than peat cores with no S. balticum (Table 1). However, the average δ¹³CH₄ emitted from the partly thawed chamber plots was −79.6 ± 0.9 ‰ in 2011, similar to or ¹³C-depleted relative to the porewater CH₄ at the same sites, indicating little CH₄ oxidation (McCalley et al., 2014). Thus, the heavier δ¹³CH₄ in the fully and partly thawed stages could also indicate increased CH₄ oxidation close to O₂-rich air at the topmost peat or in the oxidized rhizosphere at ca. 20–44 cm depth (Perryman et al., 2020; Singleton et al., 2018). In addition, during CH₄ oxidation, methanotrophs preferentially consume lighter ¹²CH₄, enriching the remaining CH₄ pool in ¹³C while producing isotopically light CO₂ (Whiticar, 1999). However, isotopically light CO₂ could also indicate less methanogenesis instead of increased methanotrophy, and rhizospheric δ¹³CH₄ enrichment can also result from faster plant ¹²CH₄ uptake and increased plant-mediated CH₄ transport (Chanton, 2005). Thus, based on the lower α_C (i.e., maintained acetoclastic methanogenesis) at 20–44 cm depth, and the strongly negative relationships between root TD and CH₄ fluxes (i.e., decreasing amount of aerenchyma is related with lower CH₄ fluxes and vice versa), it seems that the CH₄-oxidizing effect of herbaceous plants particularly in the fully thawed stage is overshadowed by increased carbon provision for methanogenesis and enhanced plant-mediated CH₄ transport.

4.3 Belowground trait relationships with CH₄ fluxes did not vary within the productive season

The relationships between the static belowground traits and CH₄ fluxes (early, middle, peak, and season median) did not vary strongly within the productive season. This may result from the lack of temporal root and rhizome sampling within the productive season. However, as an exploratory analysis, our results may indicate higher herbaceous plant-mediated CH₄ transport at the middle CH₄ flux period (Fig. 3). In contrast to our hypotheses, shrub root traits were more strongly related (i.e., higher linear regression β and R²) with peak and season median CH₄ fluxes whereas herbaceous root (TD) and rhizome (SA) traits were strongly related with CH₄ fluxes throughout the season. This could result from the lack of temporal sampling, but also from low temporal variation in coarse root diameter (diameter is used for determining root TD; see Sect. 2.5) (Picon-Cochard et al., 2012), possibly contributing to the lack of temporal trends between root diameter and CH₄ fluxes. The root TD-CH₄ flux relationships may have reflected higher CH₄-transporting herbaceous root biomass at the middle and peak CH₄ flux period (July–September; Shaver and Billings, 1977) instead of temporal TD and diameter variation. Indeed, diffusive plant-mediated CH₄ transport via Carex spp. could be high during the mid- and peak-growing season, as for example Carex rostrata (present at the partly and fully thawed stages but not captured by our aboveground measurements; see Table C1 and e.g., Hough et al., 2022) has been found to transport most CH₄ at high summer and early autumn (Ge et al., 2025; Noyce et al., 2014). Furthermore, a study conducted in Stordalen during the same productive season found that plant-mediated CH₄ transport and carbon substrate provision via graminoids was highest in August and September, further confirming these results (Mollenkopf et al., 2026). In addition, gross primary production (GPP), a proxy for plant-mediated gas transport (Knox et al., 2021), peaked in July and August (max monthly medians: 3.87 and 2.79 µmol CO₂ m⁻² s⁻¹, respectively; Fig. B1). Concurrently, CH₄ fluxes at the partly and fully thawed stages were slightly higher during nighttime (Fig. B1). Both of these GPP and CH₄ flux trends could further indicate seasonally (via higher herbaceous plant biomass) and diurnally- (via stomatal opening and closure) varying plant-mediated CH₄ transport . However, as GPP and CH₄ flux peaks did not coincide in the same hourly time periods, further research into diel CH₄ fluxes and plant-mediated CH₄ transport is needed (Fig. B1) (Knox et al., 2021). As Carex, Eriophorum and Equisetum rhizomes are perennial (Husby, 2013; Shaver, 1976), the pre-existing rhizomes may have contributed to CH₄ transport particularly in the beginning of the productive season with active root growth (esp. E. angustifolium; Shaver and Billings, 1975), resulting in strong rhizome-CH₄ flux relationships (Figs. 4, B7).

Our exploratory analyses suggested stronger shrub carbon substrate provision (SRL and diameter) at peak CH₄ flux period and across the season (season median; Fig. 3). While uncertain due to the static trait values, these results could be supported by the greater GPP at the intact and partly thawed stages (Łakomiec et al., 2021) in July and August (Fig. B1), as GPP can also be a proxy for root exudation and correlates with increased wetland CH₄ emissions (Knox et al., 2021). In addition, in a separate study within the same productive season, shrubs and graminoids were found to release more root-derived carbon at the end of the season (Mollenkopf et al., 2026). However, the biomass-weighted SRL trends were negative, which could result from species-specific SRL variation and the heterogeneous herbaceous root system (Gao et al., 2024; Picon-Cochard et al., 2012). The stronger negative trends in the beginning of the productive season may arise from increased plant carbon investment into herbaceous root growth in comparison with shrubs (Billings et al., 1978; Kummerow and Russell, 1980; Wang et al., 2016). However, within the intact and partly thawed stages and in a mixture of herbaceous and shrub plants, the SRL-CH₄ flux trends could in fact be positive (Fig. 3), but inferential statistics with larger sample sizes per thaw stage are required to investigate this further. As temporal variation in SRL and its relationship with root exudation can be negligible throughout growing season (Gao et al., 2024) and studies in (sub)arctic peatlands are lacking (Iversen et al., 2015), the negative trends between biomass-weighted SRL and CH₄ fluxes across the thaw gradient remain uncertain.

Taken together, the temporal variation in root and rhizome traits may have contributed to CH₄ fluxes but belowground trait sampling with larger sample sizes throughout the productive season is needed to investigate these relationships further. Given the uncertain herbaceous SRL and root diameter relationships, future studies could benefit from separating the roots into coarse and fine fractions, measuring the associated root exudation rates, and utilizing minirhizotrons (Iversen et al., 2012; Phillips et al., 2008) to disentangle the temporal variation in root carbon substrate provision and CH₄ transport across the permafrost thaw gradient.

5 Conclusions

Permafrost thaw alters peatland methane (CH₄) fluxes through changes in soil anoxia, nutrient and labile carbon availability, and vegetation composition. Here, we explored the relationships between plant root and rhizome traits and CH₄ fluxes along a permafrost thaw gradient in a subarctic peatland. The belowground traits reflected vegetation community shifts along the thaw gradient from shrub- to herbaceous-dominated, with decreasing shrub belowground biomass and root tissue density with increasing thaw. Root tissue density showed a strong negative correlation with CH₄ fluxes, and increasing herbaceous rhizome surface area was associated with higher CH₄ fluxes. These relationships indicated increasing plant-mediated CH₄ transport with thaw, which was supported by lower pore water CH₄ concentrations in the fully thawed stage. Simultaneously, shrub specific root length (SRL) was positively, and root diameter negatively related with CH₄ fluxes, suggesting increased labile carbon substrate provision for acetoclastic methanogenesis. The relationships did not vary strongly during the productive season, but this was likely a result of the static belowground trait data, warranting further temporal trait sampling.

We encourage further observational and experimental research particularly into the relationships between root and rhizome tissue density, SRL and CH₄ fluxes. The observed relationships could be further investigated with larger sample sizes and by combining temporal root trait sampling, root separation into fine and coarse fractions, and root exudation measurements to especially quantify the relationship between SRL and root exudation. Nonetheless, as root and rhizome trait data are scarce particularly from permafrost peatlands, and as the contribution of roots and rhizomes on plant-mediated CH₄ fluxes is understudied, this study provided valuable belowground proxies for plant-mediated CH₄ transport and fluxes. These proxies can inform future observational and experimental studies, as well as process-based wetland models by improving parameterization particularly related to plant-mediated CH₄ transport and helping fill gaps in our understanding of plant-mediated CH₄ cycling.

Appendix A: Texts

A1 Vegetation surveys

To evaluate the similarity between the CH₄ flux chamber plots and peat core plots, we conducted vegetation surveys using the point-intercept method on 19–23 July 2023. We used a 1.01 m² survey grid divided into 36 subquadrats (16.8 cm × 16.8 cm = 0.03 m²) and 25 intercept points. We lowered a dowel (60 cm long, 0.5 cm diameter) at each intercept point and noted the presence (i.e., “hit”) of each vascular plant and moss species that the dowel touched. At the peat core plots, the grid was always positioned towards northeast from the measurer and with the core in the middle of the grid. At the chamber plots, the grid was positioned so that the chamber collar was approximately in the middle of the northeast edge of the grid (due to the chamber structure, we were unable to set the collar perfectly in the middle of the grid). The ground cover (e.g., litter and water), plant and moss species were surveyed separately. For ground cover data, if the measuring dowel hit multiple classes, only the class that touched the dowel first was recorded. Species and ground cover percentages were obtained by dividing the sum of hits for a species, vegetation (herbaceous, shrub and moss) or ground cover class by the total number of intercepts (n= 25), multiplied by 100.

A2 Chamber CH₄ flux analyses

The fully thawed permafrost stage had the highest CH₄ flux (season median = 75.2 mg CH₄ m⁻² d⁻¹, IQR = 60.6 mg CH₄ m⁻² d⁻¹), followed by partly thawed stage (season median = 30.0 mg CH₄ m⁻² d⁻¹, IQR = 33.6 mg CH₄ m⁻² d⁻¹) but CH₄ fluxes at partly and fully thawed stages did not differ significantly (p>0.05). Intact stage had the lowest CH₄ fluxes (season median = 0.3 mg CH₄ m⁻² d⁻¹, IQR = 2.1 mg CH₄ m⁻² d⁻¹, p<0.01). At all thaw stages, CH₄ flux increased from the beginning to the end of the productive season: the maximum monthly median CH₄ fluxes were observed in August (max CH₄ flux: fully thawed 375.9 mg CH₄ m⁻² d⁻¹, partly thawed 160.4 mg CH₄ m⁻² d⁻¹, intact 96.4 mg CH₄ m⁻² d⁻¹), while minima were reached in June at intact (−1.9 mg CH₄ m⁻² d⁻¹) and fully thawed (10.3 mg CH₄ m⁻² d⁻¹) stages, and in May at the partly thawed stage (2.0 mg CH₄ m⁻² d⁻¹). In general, variation in CH₄ flux (coefficient of variation, CV) was larger between thaw stages (94 %; May: 109 %, August: 79 %) than within them but the intact stage had large within-stage CH₄ flux variation (148 %) which increased from May and June (125 % and 124 %, respectively) to August (165 %), while the other stages had relatively low variation (partly thawed: 37 %, fully thawed: 29 %).

Appendix B: Figures

Figure B1 Daily net ecosystem CO₂ exchange (NEE; µmol CO₂ m⁻² s⁻¹) (a) and gross primary production (GPP; µmol CO₂ m⁻² s⁻¹) (b), and hourly GPP and CH₄ fluxes during the productive season in 2023. (a) The productive season was defined as the period between the first and last passing of daily median NEE past zero (highlighted in gray; 19 May 2023 to 30 August 2023) (Körner et al., 2023). (b) Daily median GPP (based on nighttime partitioning) is shown for reference as a proxy for plant-mediated carbon substrate provision at the ecosystem scale (Knox et al., 2021). (c) Hourly GPP medians (based on daytime partitioning with light response curves for better estimation of diurnal trends). The error bars represent the ranges between 25th and 75th percentiles around the hourly medians. (d) Hourly CH₄ flux medians (in red) and means (in blue) at each permafrost thaw stage. The error bars represent the ranges between 25th and 75th percentiles around the hourly medians. Hourly means are shown in addition to medians to reflect short-term CH₄ pulses at the diurnal scale. NEE and GPP data: Lundin et al. (2023). Note that the NEE and GPP data likely mostly originate from intact and partly thawed permafrost stages and do not cover the fully thawed stage (Łakomiec et al., 2021).

Figure B2 The similarity of peat core (n= 9) and automated chamber CH₄ flux measurement (n= 9) plots based on principal component analysis (PCA). In the figure, chamber and core plots are represented by dark purple and orange points, respectively, while thaw stages are shown in different shapes (rectangle = intact permafrost, circle = partly thawed permafrost, triangle = fully thawed permafrost). The ellipses represent 95 % data ellipses for the chamber and core plots within the multivariate space defined by the first two principal components (chamber = dark purple, orange = core). The PCA was based on normalized plot herbaceous and shrub plant cover and green area (GA), mean peat temperature (mean TS) and peat or porewater pH (here, referred to as peat pH), using function PCA from package FactoMineR (Lê et al., 2008) (results shown in the figure). In addition, we used PERMANOVA with function adonis2 from package vegan (Oksanen et al., 2025) and Euclidean distances between plots within thaw stages to estimate the similarity between core and chamber plots based on the same environmental variables. Based on these assessments, the peat core plots had a larger spread in multivariate space than CH₄ flux plots mostly due to slightly higher shrub or herbaceous plant cover in one plot within each thaw stage, but the differences were not significant with the sample size n= 9 (p= 0.94). Thus, the root and rhizome measurements could be used to roughly represent belowground traits in the CH₄ flux measurement plots. However, we acknowledge that the measured root and rhizome traits do not exactly correspond to those within the CH₄ flux plots and took this into account in the result interpretation.

Figure B3 Root traits across the depth profile (0–10, 10–20, 20–30 cm) per thaw stage and plant functional type (PFT). The stacked bar colors (dark blue: herbaceous, light blue: shrub) represent the proportion of the PFT on the total biomass (a) or surface area (b), or the contribution of each PFT on the biomass-weighted trait mean (c–e). The points (light green: herbaceous, light purple: shrub) represent depth-per-plot measurements. (a–b) and (e) Root biomass, surface area and length were standardized to 10 cm depth increments. (c)–(e) Biomass-weighted trait means.

Figure B4 Root total length (km m⁻²) across thaw stages per and across plant functional types (PFT). Large points indicate thaw stage mean and small points plot-scale measurements. Colors indicate PFTs (dark grey: herbaceous, light grey: shrub, white: PFTs summed). Error bars represent 1 standard deviation of the mean. Root total length was standardized to 30 cm depth in all plots.

Figure B5 Rhizome traits across the depth profile (0–10, 10–20, 20–30 cm) per thaw stage and plant functional type (PFT). The stacked bar colors (dark blue: herbaceous, light blue: shrub) represent the proportion of the PFT on the total biomass (a) or surface area (b), or the contribution of each PFT on the biomass-weighted trait mean (c–d). The points (light green: herbaceous, light purple: shrub) represent depth-per-plot measurements. (a–b) and (e) Root biomass, surface area and length were standardized to 10 cm depth increments. (c–d) Biomass-weighted trait means.

Figure B6 Relationship between root biomass (columns a–b), surface area (columns c–d) and CH₄ fluxes (early, middle, peak and season; rows). (a, c) Relationships between root traits and CH₄ fluxes per plant functional type (PFT), where dark gray: herbaceous plants, light gray: shrubs. (b, d) Relationships between root traits and CH₄ fluxes across PFTs (biomass and surface area summed), where light blue rectangle: intact, gray circle: partly thawed, and dark purple triangle: fully thawed stage. Kendall's τ is shown for relationships where τ≥0.2 (note that due to the low sample size, p-values should be interpreted with caution). Linear regression model assumptions were not met for any of the relationships.

Figure B7 Relationship between rhizome biomass (columns a–b), tissue density (columns c–d) and CH₄ fluxes (early, middle, peak and season; rows). (a, c) Relationships between rhizome traits and CH₄ fluxes per plant functional type (PFT), where dark gray: herbaceous plants, light gray: shrubs. (b, d) Relationships between rhizome traits and CH₄ fluxes across PFTs (summed biomass, biomass-weighted tissue density), where light blue rectangle: intact, gray circle: partly thawed, and dark purple triangle: fully thawed stage. The lines represent linear regression fit used for visualization, while R² and effect size (β) are reported for linear regressions with centered and scaled CH₄ flux and trait values (when R²>0.2, model p<0.1, and linear regression model assumptions were met). Kendall's τ is shown for relationships where τ≥0.2 (note that due to the low sample size, p-values should be interpreted with caution).

Figure B8 Relationship between plant functional type (PFT)-specific plant green area (column a), plant green area summed across PFTs (column b), peat temperature (column c), peat moisture (volumetric water content VWC; column d), and CH₄ flux season values (rows). Kendall τ correlation coefficient (τ) is reported for the relationship when τ>0.2 (note that due to sample n= 9, the correlation p-values should be interpreted with caution). The different point shapes represent thaw stages (rectangle = intact permafrost, circle = partly thawed permafrost, and triangle = fully thawed permafrost). (a) Colors represent different PFTs (dark blue = herbaceous plants, light blue = shrubs).

Figure B9 The plant green area (GA; m² m⁻²) varied slightly between thaw stages. Large points represent group means, small points plot-scale measurements, and error bars ±1 standard deviation. The colors represent plant functional type (PFT) groups: dark grey is herbaceous, light grey is shrub, and white is the sum of PFT GAs. Herbaceous GA coefficient of variation (CV) was highest within intact permafrost (137 %), while shrub GA varied most within intact permafrost (109 %). Note that shrubs (Betula nana) were only found in one plot and were only 3 % of the total GA in the fully thawed permafrost stage.

Figure B10 Daily median CH₄ flux with CH₄ flux period values over the productive season at intact (a), partly thawed (b) and fully thawed (c) stages. (a–c) Red lines represent chamber-specific GAM fit (R²≥0.3), except for one chamber in a) where GAM R²<0.1, and the line represents the moving 7 d median CH₄ flux. Squares represent early, circles middle and triangles peak CH₄ flux season. See chamber GAM R² in Table C4.

Figure B11 Porewater CH₄ concentration (mM) (a), δ¹³C-CH₄ (‰) (b), CO₂ concentration (c), δ¹³C-CO₂ (‰) (d), and carbon isotope fractionation factor (α_C CO₂ $/$ CH₄) (e) from peat surface to 84 cm depth at partly and fully thawed permafrost stages. Line and point colors represent thaw stages. (c) α_C: higher values indicate larger fractionation (i.e., more hydrogenotrophy) and vice versa (i.e., more acetoclasty). Briefly, mean CH₄ concentrations increased from 0.027 ± 0.004 (standard deviation) mM to 0.533 ± 0.0007 mM at 70–74 cm depth at the partly thawed stage and from 0.017 ± 0.007 mM to 0.238 ± 0.213 mM at 80–84 cm depth at the fully thawed stage. δ¹³C-CH₄ values became enriched at ca. 30–44 cm (partly thawed ca. −56.4 ‰, fully thawed ca. −46 ‰), and this enrichment was associated with a decrease in α_C (ca. 1.048 partly thawed, ca. 1.03 fully thawed). Mean CO₂ concentrations increased from 1.758 ± 0.284 mM to 11.454±0.766 at 60–64 cm depth and 10.923 ± 1.51 at 70–74 cm depth at the partly thawed stage (with a decrease at 30–44 cm depth, with means of 3.951 and 4.430 mM), and from 1.322 ± 0.089 to 3.576 ± 1.725 mM at 80–84 cm at the fully thawed stage. δ¹³C-CO₂ became more enriched with depth particularly at the partly thawed stage (from −19.398 ‰ to −3.394 ‰).

Appendix C: Tables

Table C1 Plant cover (%) by species and herbaceous and shrub vegetation classes, litter and water surface cover (%), and peat characteristics in each chamber and peat coring (“core”) plot (as peat cores were selected to represent each individual chamber, the core and chamber plots have the same ID number). Note that the intact 1 chamber plot was on a collapsing palsa. Chamber plot WTD is the mean annual WTD across the 2007–2017 period per thaw stage (Crill et al., 2023). Values after “±” are one standard deviation of the mean. WTD = water table depth (cm; from the peat surface).

Table C2 Species included in the vegetation green area (GA) measurements per thaw stage. Mean GA (with standard deviation) is reported for each species across the three plots per thaw stage, and comprises the sum of stem and leaf green area (m² m⁻²).

Table C3 Vascular plant species growing directly on the peat cores and their leaf, stem, and flower biomass. As peat cores were selected to represent each individual chamber, the core and chamber plots have the same ID number. Note that the biomass estimates are likely underestimations due to some plant material being lost in the field. We show this table primarily for the species list. n/a: not applicable.

Table C4 Results of generalized additive models (GAM) used for determining CH₄ flux season values. Note: Intact 2 plot GAM $R_{\mathrm{adj}}^{\mathrm{2}}\mathit{<}\mathrm{0.1}$, and the CH₄ flux periods are based on moving 7 d medians. CI = 95 % confidence interval.