High peatland methane emissions following permafrost thaw: enhanced acetoclastic methanogenesis during early successional stages

Permafrost thaw in northern peatlands often leads to increased methane (CH₄) emissions, but the underlying controls responsible for increased emissions and the duration for which they persist have yet to be fully elucidated. We assessed how shifting environmental conditions affect microbial communities and the magnitude and stable isotopic signature (δ¹³C) of CH₄ emissions along a thermokarst bog transect in boreal western Canada. Thermokarst bogs develop following permafrost thaw when dry, elevated peat plateaus collapse and become saturated and dominated by Sphagnum mosses. We differentiated between a young and a mature thermokarst bog stage (∼ 30 and ∼ 200 years since thaw, respectively). The young bog located along the thermokarst edge was wetter, warmer, and dominated by hydrophilic vegetation compared to the mature bog. Using high-throughput 16S rRNA gene sequencing, we show that microbial communities were distinct near the surface and converged with depth, but fewer differences remained down to the lowest depth (160 cm). Microbial community analysis and δ¹³C data from CH₄ surface emissions and dissolved gas depth profiles show that hydrogenotrophic methanogenesis was the dominant pathway at both sites. However, mean δ¹³C-CH₄ signatures of both dissolved gas profiles and surface CH₄ emissions were found to be isotopically heavier in the young bog (−63 ‰ and −65 ‰, respectively) compared to the mature bog (−69 ‰ and −75 ‰, respectively), suggesting that acetoclastic methanogenesis was relatively more enhanced throughout the young bog peat profile. Furthermore, mean young bog CH₄ emissions of 82 mg CH₄ m⁻² d⁻¹ were ∼ 3 times greater than the 32 mg CH₄ m⁻² d⁻¹ observed in the mature bog. Our study suggests that interactions between the methanogenic community, hydrophilic vegetation, warmer temperatures, and saturated surface conditions enhance CH₄ emissions in young thermokarst bogs but that these favourable conditions only persist for the initial decades after permafrost thaw.

1 Introduction

Methane (CH₄) emissions in northern peatlands are typically thought to be driven by environmental and ecological conditions such as temperature, water table position, and vegetation community (Bellisario et al., 1999). However, CH₄ emissions are ultimately the result of microbial activity, and understanding the interactions between environmental conditions and microbial processes is key to understanding the impact of disturbances on peatland CH₄ emissions. Increased disturbances such as permafrost thaw are transforming northern-latitude peatlands (Helbig et al., 2016) through the disruption of the frozen landscape and environmental conditions responsible for the regional accumulation of large peatland carbon (C) stores. Rapidly rising northern air temperatures (Mudryk et al., 2018) are predicted to lead to widespread gradual thawing of permafrost (Schaefer et al., 2011) and subsequent thermokarst development in high-C-density permafrost peatlands (Olefeldt et al., 2016). Thermokarst formation in ice-rich permafrost peatlands is characterized by ground subsidence and surface inundation (Camill, 1999). This exposes previously frozen C to anaerobic microbial decomposition and potential mineralization into greenhouse gases (Schuur et al., 2015). Redox conditions following thermokarst formation are an important control of decomposition, with 3–4 times greater C mineralization occurring as aerobic respiration compared to anaerobic respiration (Schädel et al., 2016). Increased emissions of methane (CH₄) due to thermokarst formation are projected to result in a positive feedback with climate warming (Turetsky et al., 2020). However, the magnitude of peatland CH₄ emissions and the metabolic pathways responsible for these emissions in response to permafrost thaw remain uncertain, as does the period for which these conditions and emissions persist.

Methanogenesis, conducted by methanogenic archaea belonging to phylum Euryarchaeota, is one of the most prominent microbial processes contributing to the anaerobic decomposition of organic matter in water-logged permafrost soils (Cai et al., 2016; Knoblauch et al., 2018). Methanogenesis occurs primarily via two pathways: acetoclastic methanogenesis and hydrogenotrophic methanogenesis (Whiticar et al., 1986; Whiticar, 1999). Acetoclastic methanogenesis involves the cleavage of acetate into CH₄ and CO₂ and, when considering these two species, causes less apparent fractionation than the hydrogenotrophic methanogenesis pathway. This results in acetoclastic methanogenesis yielding comparatively isotopically heavy δ¹³C-CH₄ (δ¹³C = −65 ‰ to −50 ‰). The reduction of CO₂ and H₂ in hydrogenotrophic methanogenesis typically produces CH₄ lighter in ¹³C (δ¹³C = −110 ‰ to −60 ‰) (Hornibrook et al., 1997, 2000). While the two pathways are stoichiometrically equal (Conrad, 1999; Corbett et al., 2013), the activity of acetoclastic and hydrogenotrophic methanogens is governed by different extrinsic controls (Bridgham et al., 2013).

Hydrogenotrophic methanogenesis is thought to be the main pathway of CH₄ formation in northern peatlands (Hornibrook et al., 1997; Galand et al., 2005). However, the acetoclastic pathway can dominate in the upper layers of more minerotrophic, nutrient-rich peatlands (Popp et al., 1999; Chasar et al., 2000) where there are sufficient levels of acetate (Ye et al., 2012). During the initial decades following thaw, surface runoff of nutrients from surrounding intact peat plateaus (Keuper et al., 2012, 2017) and increased connectivity to regional hydrology (Connon et al., 2014) can result in more minerotrophic conditions. Such shifts in hydrology, temperature, nutrients, redox conditions, and vegetation communities following permafrost thaw have been shown to increase the prevalence of acetoclastic methanogenesis and CH₄ emissions (Hodgkins et al., 2014; McCalley et al., 2014). However, this potential post-thaw enhancement of acetoclastic methanogenesis needs to be considered in the context of the existing methanogenic community that developed in the peat profile before thaw. For example, historical environmental conditions have been shown to have a legacy effect on the methanogenic community following thaw and can therefore be a key constraint on methanogenic community structure and activity post-thaw (Holm et al., 2020; Lee et al., 2012). Overall, an understanding of the methanogenic community's response following thaw to shifts in both surface conditions and exposure to previously frozen organic matter is key to estimating CH₄ emissions from thermokarst peatlands.

Environmental conditions following permafrost thaw in peatlands are characterized by a drastic shift in water table position and increased wetness, increased soil temperatures, and a change in vegetation community associated with increased labile inputs (Beilman, 2001; Burd et al., 2020; Camill, 1999). These shifts may provide optimal conditions for CH₄ production and emissions, particularly in the initial decades following thaw. Peatland CH₄ emissions are constrained by water table position (Huang et al., 2021; Strack et al., 2004), and surface inundation leads to increased CH₄ emissions (Tuittila et al., 2000). Methane production and emissions are positively influenced by soil temperatures (Hopple et al., 2020; Olefeldt et al., 2017), and peatland CH₄ emissions have been shown to increase when both water table position and temperatures are high (Grant, 2015). The colonization of vegetation associated with fresh, labile inputs has also been shown to increase both the magnitude and temperature sensitivity of CH₄ emissions in peatlands (Leroy et al., 2017; McNicol et al., 2019). As such, many studies have focussed on the relationship between water table position, soil temperature, and vegetation communities in determining CH₄ fluxes following thaw (Johnston et al., 2014; Turetsky et al., 2007; Wickland et al., 2006). However, while these environmental conditions are key drivers of CH₄ emissions, they are unable to fully account for the variability in permafrost peatland CH₄ emissions (Juottonen et al., 2021; Kuhn et al., 2021). Some of this unaccounted variance may be in part explained by microbial activity, as changes in the composition and abundance of methanogenic community members can contribute significantly towards peatland CH₄ emissions (Fritze et al., 2021). Relatively few studies have assessed how shifts in environmental conditions and ensuing changes in methanogenic community structure influence CH₄ emissions following thaw (McCalley et al., 2014), an interaction that may be significant at both the local and circumpolar scale.

In this study, we assess the impact of permafrost thaw on peatland methanogenic community composition and CH₄ emissions along a space-for-time thaw gradient that includes an intact peat plateau and an adjacent thermokarst bog with areas that have thawed ∼ 30 and ∼ 200 years ago (herein referred to as young bog and mature bog, respectively). Thermokarst formation has resulted in distinct environmental conditions at each stage along this thaw gradient. We herein define these distinct environmental conditions as water table position and surface wetness, soil temperatures, and vegetation community. Along this gradient, we assessed methanogenic community structure down to 160 cm. We hypothesize that (1) shifting environmental conditions along the permafrost thaw gradient result in a successional microbial community and a restructuring of the methanogenic community, and (2) the warmer conditions and hydrophilic vegetation community in the young bog, along with the exposure of previously frozen peat, will result in a greater relative abundance of acetoclastic methanogens throughout the depth profile and subsequently greater overall CH₄ emissions. In the young bog and mature bog, we measured the concentration and δ¹³C signature of dissolved CH₄ and CO₂ down to 245 cm, as well as the rates and δ¹³C signature of both CH₄ and CO₂ land–atmosphere fluxes. The combined approach of measuring dissolved gas depth profiles and surface emissions, in tandem with assessing the structure of the methanogenic community along a depth profile, allows us to determine how changing environmental conditions following thaw impacts methanogenic pathways and community composition. Utilizing this approach, we can subsequently gain further insight into how long elevated surface CH₄ emissions may persist post-thaw. Furthermore, this approach highlights that while environmental conditions are important in determining CH₄ emissions, microbial community composition and changes in the methanogenic community structure are likely to significantly influence CH₄ emissions following thaw.

2 Methods

2.1 Study site and design

The Lutose peatland study site (59.5^∘ N, 117.2^∘ W; Fig. 1) is located on the Interior Plains of western Canada within the zone of discontinuous permafrost (Brown et al., 1997; Heginbottom et al., 1995). The climate is continental with a monthly average summer high temperature of 16.1 ^∘C (July), winter low of −22.8 ^∘C (January), and annual average air temperature of −1.8 ^∘C (Climate-Data.org, 2019 – data from site located ∼ 50 km south of Lutose). Annual average precipitation is 391 mm, of which three-quarters fall as rain between May and September. In the discontinuous permafrost zone of the Interior Plains in boreal western Canada, ∼ 40 % of the landscape is covered by permafrost peatlands that have between 2 and 6 m deep peat deposits (Gibson et al., 2018; Vitt et al., 2000). The peatland complexes in this area are a fine-scale mosaic of permafrost peat plateaus, as well as permafrost-free ponds, fens, and bogs (Zoltai, 1993; Bauer et al., 2003; Vitt et al., 2000; Pelletier et al., 2017), and they are similar to those found in the Hudson Bay Lowlands (Kuhry, 2008) and Alaska (Jones et al., 2017). The Lutose peatland complex is representative of the peatlands found in the discontinuous permafrost zone of the Interior Plains in western Canada (Heffernan et al., 2020). The site has 5–6 m deep peat and has transitioned through multiple developmental stages since it began accumulating organic matter ∼ 8800 years ago. It transitioned from a marsh, through a fen and a bog stage prior to permafrost aggradation ∼ 1800 years ago (Heffernan et al., 2020). Peatlands in the Interior Plains in western Canada are one of the three largest stores of organic carbon found in peatlands within the permafrost zone, the other two being the Hudson Bay Lowlands and the West Siberian Lowlands (Hugelius et al., 2020; Olefeldt et al., 2021). Within the sporadic and discontinuous permafrost zone of our study region >15 % of the total peat plateau area has thawed and formed thermokarst bogs in the last 30 years (Baltzer et al., 2014; Gibson et al., 2018). Projections for this area suggests total permafrost lost from plateaus by 2050 (Chasmer and Hopkins, 2017).

Figure 1 Lutose peatland site location and study design. (a) Site location (Lutose, Alberta, Canada; 59.5^∘ N, 117.2^∘ W) in boreal western Canada. Green shading represents permafrost zonation (Brown et al., 1997), and brown shading represents areas with >25 % permafrost peatland (histels) extent (Hugelius et al., 2014). (b) GeoEye satellite image of study site (image from https://zoom.earth/, last access: 20 November 2018) at 0.46 m resolution. The circle represents the area where sampling took place. (c) Aerial image of study transect, locations of peat and dissolved gas sampling in the plateau (blue diamond), young bog (green diamond), and mature bog (orange diamond), and area where collars for gas flux measurements were located in the young bog (YB, green) and mature bog (MB, orange) (aerial photo credit: Olefeldt, David). (d, e) Surface vegetation in the mature bog and young bog and (f) soil profile of thaw transect based on Heffernan et al. (2020). The transition to post-thaw peat occurs at 29 and 71 cm in the young bog and mature bog, respectively. Peat (core) and pore water (pore water peepers), including microbial community, sampling depth profile for 0–160 cm shown as white column with diagonal black lines. Dissolved gas (diffusive samplers) sampling depth profile for 0–245 cm shown as light green column. Soil temperature depth profile for 0–250 cm shown as purple column. Average water table depth shown as dashed blue line.

The studied transect represents a space-for-time gradient of permafrost thaw that includes three thaw stages: a permafrost peat plateau and a young (∼ 30 years since thaw) and a mature (∼ 200 years since thaw) part of an adjacent thermokarst bog. The timing of permafrost thaw was previously determined by ¹⁴C-dating the shift in macrofossil vegetation indicative of thaw, at 29 cm in the young bog and at 71 cm in the mature bog (Fig. 1f) (Heffernan et al., 2020). The peat plateau has an active layer thickness of ∼ 70 cm, and its surface is raised 1–2 m above the adjacent thermokarst bog due to the presence of excess ground ice, resulting in relatively dry surface conditions in which the water table generally follows the deepening of the seasonally thawed peat layer (Zoltai, 1972). This thaw stage is characterized by a stunted, open black spruce (Picea mariana) canopy and ground cover of lichens (Cladonia spp.), Sphagnum fuscum hummocks, and low-lying ericaceous shrubs as is characteristic of the peat plateaus in the area (Vitt et al., 1994). The young bog stage is narrow (<5–10 m wide) and is located next to the actively thawing area of the peat plateau. The young bog has an average growing season water table position of 1.3 ± 4.9 cm below the peat surface. These inundated conditions result in the dominance of a hydrophilic vegetation community (Fig. 1e) consisting of Sphagnum riparium, bog-sedge (Carex limosa), and rannoch rush (Scheuchzeria palustris). The mature bog is ∼ 10–15 m from the young bog and is drier compared to the young bog, with an average growing season water table position of 22.9 ± 9.3 cm below the surface. The dominant vegetation reflects these drier conditions and consists of Sphagnum fuscum, Sphagnum magellanicum, leather leaf (Chamaedaphne calyculata), cloudberry (Rubus chamaemorus), Eriophorum vaginatum tussocks, and some black spruce (Picea mariana) regrowth (Fig. 1d). The mature bog is located >10–20 m from the thawing plateau edge.

2.2 Site preparation and monitoring of environmental conditions

The Lutose peatland study site was established in 2015, and a boardwalk was constructed to minimize disturbances along the peat plateau–thermokarst bog transect. Three collars for surface greenhouse gas flux (39 cm diameter) measurements were permanently installed to a depth of 20 cm in both the young and mature thermokarst bog stages. The top of each collar was aligned with the peat surface. PVC wells (2 cm diameter) were installed directly next to each collar and were used to manually monitor the water table position during each gas flux measurement. We monitored soil temperature (^∘C) at 10, 30, 50, 75, 100, 150, 200, and 250 cm every 30 min from May to September 2018 using permanently installed loggers (Onset Computer Hobo 8k Pendant, Bourne, MA, USA) in both thermokarst bog stages. Temperature depth profiles were established centrally among collars in each thermokarst bog stage and in areas that had similar vegetation, water table position, and distance from the thawing edge as the collars.

Custom made plexiglass pore water suction (Heffernan et al., 2021) and diffusive equilibration gas sampling devices (Knorr et al., 2009) were installed in July 2016 in the young and mature bogs. These devices were installed in both thermokarst bog stages ∼ 1 m from the nearest flux measurement collar. Pore water suction devices were installed to a depth of 160 cm and consisted of 15 sampling depths, with each sampling depth connected to the surface via silicone tubing. This allowed for repeated non-destructive pore water sampling. Three diffusive gas sampling devices were installed in each thermokarst bog stage, where two collected dissolved soil gas samples from 5 to 95 cm deep and a third from 115 to 245 cm. Each diffusive gas sampler consisted of a PVC pipe with a 10 cm long sampling section centred at each sampling depth. Sampling sections consisted of ∼ 2 m of silicon tubing (3 mm i.d., 5 mm o.d.) wrapped around the PVC pipe and kept in place by PVC spacers at the top and bottom of each interval. Silicone tubes were sealed at one end, whereas the other end was connected to polyurethane tubing (1.8 mm i.d.) that ran back up inside the PVC tube to reach the peat surface where it was sealed with a three-way stopcock. Silicone tubing has been shown to be permeable to gases such as CO₂ and CH₄ within a number of hours while remaining impermeable to water, making it suitable for sampling of dissolved soil gases (Kammann et al., 2001).

2.3 Pore water chemistry and peat enzyme activity

Pore water dissolved organic matter (DOM) chemistry and peat enzyme activity presented in this study have previously been published (Heffernan et al., 2021) and are briefly described here. Pore water samples for DOM chemistry were taken monthly from May to September 2018 using the previously described pore water suction devices in the young bog and mature bog. Three 60 mL samples were taken from all 15 measurement depths by applying a vacuum at the surface and collecting water with syringes via a three-way stopcock. Each water sample was immediately filtered through 0.7 µm pore size glass fiber filters (GF/F Whatman) into two acid-washed amber glass bottles, with one sample acidified with 0.6 mL 2N HCl to prevent further microbial activity. Pore water samples were transported in a cooled container and stored at 4 ^∘C prior to analysis. Pore water DOM was analyzed for pH, phosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-}$; µg L⁻¹), dissolved organic carbon (DOC; mg L⁻¹), and total dissolved nitrogen (TDN; mg L⁻¹) concentrations, phenolic contents, specific UV absorbance at 254 nm (SUVA; L mg C⁻¹ m⁻¹; Weishaar et al., 2003), and spectral slope between 250 and 465 nm (S_250−465, nm⁻¹; Helms et al., 2008). SUVA and S_250−465 values are used to indicate aromaticity, with high SUVA indicating a high aromatic content and lower S_250−465 indicating low molecular weight and decreasing aromaticity (Hansen et al., 2016).

Peat cores extracted to a depth of 160 cm were stored at 4 ^∘C for less than 1 week in the laboratory before homogenization to determine potential soil enzyme activities. We performed hydrolytic enzyme assays for four enzymes – phosphatase, β-N-glucosaminidase, β-glucosidase, and β-cellobiosidase – using fluorogenic 4-methylumbelliferone labelled substrates (Dunn et al., 2014). We assayed oxidative enzyme activity by measuring laccase activity using syringaldazine (Criquet et al., 2000; Jassey et al., 2012). We summarized the activity of all enzymes using a multi-functionality index based on z scores (Allan et al., 2015; Heffernan et al., 2021).

2.4 Surface land–atmosphere gas fluxes

We measured surface land–atmosphere greenhouse gas fluxes (CH₄ and carbon dioxide; CO₂) monthly from May to September 2018 at the three collars in each peatland stage using the static chamber method (Carroll and Crill, 1997). The chamber used to capture land–atmosphere fluxes was a transparent cylindrical Plexiglass chamber with a basal area of 0.12 m², height of 0.40 m, and volume of 47.8 L. The chamber was equipped with three fans (Micronel D341T012GK-2 ventilator, BEDEK GmbH, Dinkelsbühl, Germany) to mix air during measurements and a temperature sensor (Hobo RH Smart Sensor, S-THB-M002, Onset Computer, Bourne, USA) that was shaded from direct sunlight (Burger et al., 2016). An airtight seal was formed between the chamber and collar by pouring water in a ∼ 1.5 cm deep well around the upper circumference of each collar. Land–atmosphere fluxes of CO₂ (ecosystem respiration) and CH₄ were captured simultaneously in darkened conditions by covering the chamber with a reflective shroud. Gas concentrations were determined at a temporal resolution of 1 s using an Ultraportable Greenhouse Gas Analyzer (Los Gatos Research, CA, USA), and real-time fluxes were monitored using the VNV^® Viewer (RealVNC^® Limited, UK) application with an iPad mini 2 (Apple Inc.).

The rates of CH₄ and CO₂ land–atmosphere fluxes (Flux) were calculated using the following ideal gas law:

$$\begin{array}{}\text{(1)}& \mathit{\text{Flux}}=\mathit{\text{slope}}{\displaystyle \frac{P\cdot V}{R\cdot T\cdot A}},\end{array}$$

where slope is the linear rate of change in gas concentration (µmol mol⁻¹ s⁻¹) over the measurement period inside the chamber, P is an atmospheric pressure (atm) constant of 0.96 atm, V is chamber volume (L), R is the universal gas constant (L atm K⁻¹ mol⁻¹), T is the average temperature (K) inside the chamber during the measurement, and A is the chamber basal area (m²). Chamber closure for each flux measurement was 5 min with the first 2 min discarded to ensure fluxes (i.e., change in concentration over time) with R²>0.75. We report CO₂ fluxes in grams CO₂ per square meter per day (g CO₂ m⁻² d⁻¹) and CH₄ fluxes in milligrams CH₄ per square meter per day (mg CH₄ m⁻² d⁻¹), with positive values indicating fluxes to the atmosphere. To quantify the proportion of C being emitted as CH₄, we standardized our CO₂ and CH₄ fluxes per gram C emitted. The proportion of C emitted as CH₄ (CH₄ : C emissions) was calculated as

$$\begin{array}{}\text{(2)}& {\mathrm{CH}}_{\mathrm{4}}:\mathrm{C}\mathrm{emissions}={\displaystyle \frac{{\mathrm{CH}}_{\mathrm{4}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}}{{\mathrm{CH}}_{\mathrm{4}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}+{\mathrm{CO}}_{\mathrm{2}}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}}}.\end{array}$$

2.5 δ¹³C signature of CH₄ emissions

We assessed the δ¹³C-CO₂ and δ¹³C-CH₄ signatures of ecosystem respiration (CO₂) and CH₄ emissions. This was done similarly to regular measurements of CO₂ and CH₄ fluxes but using a smaller, opaque chamber of 31.1 L and discrete syringe samples for δ¹³C analysis in combination with the continuous monitoring of gas concentrations described above. Gas syringe samples were taken using a 20 mL syringe via a three-way stopcock placed between the sealed chamber and gas inlet port on the Ultraportable Greenhouse Gas Analyzer. Gas samples were then injected into a 37.5 mL sealed glass vial that had been flushed with nitrogen gas prior to sealing. Chamber enclosure time ranged from 30 to 50 min with four to five samples being taken during this time. Samples were taken either every 10 min or once a minimum change in CO₂ (30 µmol mol⁻¹) and CH₄ (1 µmol mol⁻¹) concentrations was observed. An atmospheric gas sample was used as a time-zero measurement when assessing the change in concentration over time. Glass vials containing samples were stored at 4 ^∘C until analysis. These measurements were taken in September and October 2016 from one collar in both the young and mature bogs, with each collar measured twice.

We measured the δ¹³C values of gas samples from both the chamber fluxes and atmospheric background. To assess whether the gas concentration of each sample fit within the measurement range required for δ¹³C analysis, we measured CO₂ and CH₄ concentrations using 1–3 mL from each vial. Following these concentration measurements, the remaining sample (17–19 mL) was diluted with nitrogen gas to a final volume of 20 mL and injected into a small sample introduction module (SSIM; Picarro, California, USA) system to measure δ¹³C signatures. The δ¹³C-CO₂ and δ¹³C-CH₄ signature was measured in-line with a cavity ring-down spectrometer (G2201-L, Picarro, California, USA) that had been calibrated using certified standards.

We then used the time series of δ¹³C-CH₄ and CH₄ concentrations to estimate the δ¹³C-CH₄ signature of the CH₄ released to the atmosphere using Keeling plots (Keeling, 1958). Using this approach, the δ¹³C-CH₄ signature of gas in each sample is plotted on the y axis against the inverse of CH₄ gas concentrations ($\mathrm{1}/\left[{\mathrm{CH}}_{\mathrm{4}}\right]$). The y-axis intercept of the linear regression represents the mean isotopic signature of the CH₄ source (Fisher et al., 2017). While fractionation during diffusive transport may influence these estimates, it has been shown in similar systems to be of minor importance compared to other contributing processes (Preuss et al., 2013; Nielsen et al., 2019).

2.6 Dissolved gas depth profiles

Dissolved gas samples were collected using diffusive equilibration gas sampling devices. Samples were taken from the following 15 depths: every 10 cm down to 95 cm starting at 5–15 cm and then at 115, 140, 165, 195, and 245 cm. Once a month from May to September 2018 a ∼ 7 mL gas sample was drawn from each depth using a 10 mL plastic syringe. These gas samples were immediately injected into a 10 mL sealed glass vial that had been flushed with nitrogen gas prior to sealing and then were stored at 4 ^∘C until analysis. A total of 214 CO₂ and 211 CH₄ dissolved gas concentration measurements were made by injecting 1–3 mL of gas into a gas chromatograph with an FID (flame ionization detector) and CO₂ methanizer (8610C Gas Chromatograph, SRI Instruments, California, USA). We measured δ¹³C-CO₂ and δ¹³C-CH₄ signatures using the previously mentioned cavity ring-down spectrometer and SSIM system. As with surface chamber gas samples, dissolved gas samples were diluted with N₂ to 20 mL. However, dissolved gas concentrations were considerably higher than gas concentrations found in the surface chambers, and some were well above the optimal concentration range required for accurate δ¹³C analysis for the SSIM system even after dilution. To fit within measurement range of the system, further dilution resulted in CO₂ concentrations below detectable limits. As such, we were able to obtain 90 and 75 measurements of δ¹³C-CH₄ in the young and mature bogs, respectively, and 93 measurements of δ¹³C-CO₂ in both.

We used the δ¹³C-CO₂ and δ¹³C-CH₄ signature of each gas sample to calculate the apparent fraction factor α_c, where α_c= [¹³C-CO₂+ 1000] $/$ [¹³C-CH₄+1000]. The α_c can serve as an isotopic indicator of the pathway of methanogenesis, with typical values of 1.060–1.090 observed for hydrogenotrophic methanogenesis and 1.040–1.060 for acetoclastic methanogenesis (Chanton et al., 2005).

2.7 Peat and pore water sample collection for microbial community composition analyses

Microbial community composition was characterized in both peat and peat pore water samples from depths between 0 and 160 cm in the young bog and mature bog. Focusing on peat samples, microbial community composition in the active layer of the peat plateau was assessed from depths between 0 and 30 cm. Peat cores were extracted in June and September 2018. Near-surface cores were extracted using a cutting tool down to 30 cm depth in the peat plateau and young bog and 50 cm deep in the mature bog. Surface cores were limited to 30 cm in the plateau due to the presence of ground ice during sampling in June. Surface core depths differed between the young bog and mature bog due to differences in the water table position. Deeper core sections (down to 160 cm) in the young bog and mature bog were extracted using a Russian peat corer (4.5 cm i.d., Eijkelkamp, Giesbeek, The Netherlands). Cores were extracted from two boreholes located ∼ 20 cm apart, alternating between boreholes to avoid disturbance contamination from the 10 cm corer tip during the coring process. To do so, 50 cm long core sections were taken alternatively from each borehole, with each core having a 10 cm overlap with the previous core taken from the adjacent borehole. In the field, immediately after the entire core was extracted, cores were divided into 15 subsections. The first two subsections contained peat from 0 to 5 cm and 5 to 10 cm, followed by 10 cm increments down to 120 cm, and two further subsections from 130 to 140 and 150 to 160 cm. Peat from each interval was sub-sampled using sterilized forceps, placed directly into Whirl-Pak^® bags, and frozen within 3 h of sampling for transportation back to the laboratory. Once samples reached the laboratory, they were frozen at −80 ^∘C until analysis.

We also sampled peat pore water at all 15 peat sampling depths in September 2018 from the pre-installed pore water suction sampling devices mentioned above. We extracted 60 mL pore water samples by applying a vacuum at the surface and collecting water with new plastic 60 mL syringes. Pore water was immediately filtered through sterile 0.2 µM pore size, polyvinylidene difluoride (PVDF) membrane Sterivex filters (MilliporeSigma). Microbial cells were retained on the filter, and remaining porewater in the Sterivex was removed via extrusion using a 60 mL sterile syringe. Sterivex filters were then immediately flash-frozen at −80 ^∘C in a liquid nitrogen dry shipper to preserve microbial community members until analysis could take place.

2.8 DNA extraction

Genomic DNA was extracted from all peat and pore water samples using the DNeasy PowerSoil kit (Qiagen) and the PowerWater DNeasy kit (Qiagen), respectively, to assess the differences in microbial community structure. Extraction of DNA from both sample types was followed as described by the manufacturer (Qiagen), with two modifications: (i) for peat samples, prior to mechanical lysis using bead beating, the prepared samples were chemically lysed by incubation at 70 ^∘C for 10 min in the provided lysis solution, and (ii) Sterivex (pore water) samples were incubated with rotation at 37 ^∘C following addition of lysis buffer. These modifications were made to increase total DNA yield. The amount of isolated DNA from each sample was then determined using a Qubit fluorometer (model 2.0, using the 1× HS dsDNA kit), with concentrations ranging between ∼ 0.1 and 22.4 ng µL⁻¹. This extracted DNA served as the template for polymerase chain reaction (PCR) analyses described below.

2.9 Sequencing and computational analyses

We amplified 16S rRNA genes using universal prokaryotic primers 515F (Parada et al., 2016) and 926R (Quince et al., 2011). Each primer also contained a six-base index sequence for sample multiplexing (Bartram et al., 2011). The PCR mix (25 µL total volume) contained 1 × Q₅ reaction buffer, 0.5 µM forward primer, 0.5 µM reverse primer, 200 µM dNTPs, 0.500 U Q₅ polymerase (New England Biolabs, Ipswich, MA, USA), and 2.5 µL of genomic template. Genomic extracts with DNA concentrations greater than 2 ng µL⁻¹ were diluted 1:100 in nuclease-free water. The PCR was performed as follows: 95 ^∘C for 3 min, 35 cycles of 95 ^∘C for 30 s, 60 ^∘C for 30 s, 70 ^∘C for 1 min, and a final extension of 70 ^∘C for 10 min. Pooled 16S rRNA gene amplicons were purified using NucleoMag beads, and a 4.5 pM library containing 50 % PhiX Control v3 (Illumina, Canada Inc., NB, Canada) was sequenced on a MiSeq instrument (Illumina Inc., CA, USA) using a 2×250 cycle MiSeq Reagent Kit v3 (Illumina Canada Inc) at the Molecular Biology Service Unit (MBSU, University of Alberta). The MiSeq reads were demultiplexed using MiSeq Reporter software version 2.5.0.5. Each read pair was assembled using the paired-end assembler for Illumina sequences (PANDAseq; Masella et al., 2012) with a quality threshold of 0.9, dictating that 90 % of overlapping reverse and forward reads must match in order to assemble reads into read pairs. Assembled reads were analyzed using the Quantitative Insights Into Microbial Ecology II pipeline (QIIME2; Boylen et al., 2020). Sequences were clustered into amplicon sequence variants (ASVs) with chimeric sequences, singletons, and low abundance ASVs removed using DADA2 (Callahan et al., 2019). All representative sequences were classified with the Greengenes reference database, using the most recent release (version 13.8; McDonald et al., 2012). Although Greengenes is not updated as frequently as the SILVA database, we chose to use it to classify our ASVs as a comparison of both databases revealed that they captured a similar number of archaea (total of 51 187 methanogenic read counts attributed to SILVA versus 51 141 methanogenic read counts attributed to Greengenes). The taxonomic resolution between both databases was also similar, identifying the same kinds of phyla, families, and genera, as well as methanogens (e.g., Methanoregula, Methanosarcinales, etc.). Given these similarities, and the fact that methanogen nomenclature has not changed significantly over time, we ultimately chose to use Greengenes because it was able to resolve more methanogenic families belonging to Methanocellales and Methanomassiliicoccaceae particularly compared to SILVA. The Greengenes database is also still commonly used to explore methanogenic archaeal communities in current literature (Vanwonterghem et al., 2016; Lin et al., 2017; Carson et al., 2019). Furthermore, since 1021 methanogenic reads were captured per sample, on average, using Greengenes and are comparable to other studies (Vishnivetskaya et al., 2018; Holm, et al., 2020), we believe that our approach is sufficient for covering methanogen diversity.

2.10 Statistical analyses

All statistical analyses were carried out in R (Version 3.4.4, R Core Team, 2015) using the nlme, vegan, factoextra, ggplot2, VariancePartition, and ggpubr packages (Pinheiro et al., 2017; Oksanen et al., 2013; Kassambara and Mundt, 2017; Wickham, 2009; Hoffman and Schadt, 2016; Kassambara, 2018). For analysis of variance (ANOVA), distribution of the data was inspected visually for normality along with the Shapiro–Wilk test. We tested homogeneity of variances using the car package and Levene's test (Fox and Weisberg, 2011). We report uncertainty as ± 1 standard deviation, except for land–atmosphere greenhouse gas fluxes which we report as ± 95 % confidence intervals. We here define the statistical significance level at 5 %.

We used ANOVAs and Bonferroni post hoc tests on linear mixed-effects models to address our second hypothesis and to evaluate significant differences and seasonal trends in greenhouse gas fluxes and dissolved gas depth profiles. We performed these tests to assess whether thaw stage (young bog or mature bog) influenced greenhouse gas fluxes and dissolved gas depth profiles. This approach was used to test for significant differences in CH₄ fluxes, ratio of CH₄ : C emissions, and source ¹³C-CH₄ signature intercepts of Keeling plots between young bog and mature bog stages. In each linear mixed-effects model, sampling month and peatland stage were defined as fixed effects, whereas sampling collar was defined as a random effect. Similarly, we tested for significant differences between the young and mature bog depth profiles with respect to dissolved CH₄ and CO₂ concentrations, δ¹³C-CH₄ and δ¹³C-CO₂ values, α_c values, and pore water chemistry. In these models, sampling month and peatland stage were defined as fixed effects, while sample depth was defined as a random effect.

Following microbial 16S rRNA gene sequencing, sample reads were rarefied to the lowest read count of 28 129 for all subsequent analyses. These sequences represent whole microbial community data that were used to determine whether there was evidence of changes in microbial community structure representing the successional peatland stages following permafrost thaw throughout the 160 cm depth peat profile. In addition, to address our first hypothesis, we assessed differences in community composition across both peat and pore water to determine whether seasonality impacted microbial community structure in both sample matrices. Here, Bray Curtis dissimilarity matrices for overall microbial community data were used, at 999 permutations, to identify distinct groupings assessed at the 95 % confidence interval in NMDS (non-metric multidimensional scaling) ordinations. These distinct groupings were further evaluated for significance using the non-parametric permutational analysis of variance (PERMANOVA) test.

To further test our first hypothesis, methanogens were selected at the order level from our whole community data using Greengenes-assigned taxonomy. Utilizing their assigned taxonomy, the pathways through which identified methanogens conduct methanogenesis were determined by comparing our findings with the literature (Berghuis et al., 2019; Stams et al, 2019; Kendall and Boone, 2006; Zhang et al., 2020). Focusing on the methanogenic community allowed us to specifically assess how permafrost thaw affects the microbial community responsible for CH₄ production and net CH₄ emissions following thaw. We utilized our methanogenic community data to construct redundancy analyses (RDAs) and relative abundance bar plots. RDAs were conducted using a Hellinger-transformed methanogenic community. Explanatory variables (i.e., dissolved concentrations of CO₂, CH₄, DOC, temperature, enzymatic activity estimate, thaw stage, depth, and distance to water table) were scaled about the mean. These explanatory variables had variance standardized, were checked for collinearity (parameters with variance inflation value >10 were removed), and were selected for significance using backward selection, set at 1000 permutations. The significance of the RDA model and of each axis was tested using ANOVAs, set at 999 permutations. Variance partitioning analyses were conducted to assess the contribution of significant environmental parameters (i.e., thaw stage and distance to water table) on the structuring of the Hellinger-transformed methanogenic community. Distance from water table reflects the distance (in cm) a certain sample is from the water table in different stages of thaw (young bog and mature bog). Due to the smaller size of our methanogenic community relative to the total community and the lack of some data at certain depths, we combined pore water and peat samples together for these analyses. Relative abundance, which measures how common or rare a particular microorganism is relative to the entire microbial community, of methanogenic orders related to acetoclastic or hydrogenotrophic methanogenesis processes was plotted according to depth. Significant differences in methanogenic community composition between depths were assessed using the non-parametric Kruskal–Wallis test with a Benjamini–Hochberg correction for multiple comparisons after running a Wilcox rank sum test.

3 Results

3.1 Site environmental conditions

The young bog was wetter and warmer than the mature bog throughout the May–September 2018 study period. In June, following snowmelt, the water table was at its highest at 2.2 ± 0.6 cm above the surface in the young bog. The highest water table position in the mature bog was 17.5 ± 1.9 cm below the peat surface and observed in July. The water table dropped during the season and in September was 5.7 ± 2.2 and 27.3 ± 1.2 cm below the peat surface in the young bog and mature bog, respectively. In the plateau, the seasonally thawed layer gradually deepened during the growing season, with an active layer depth of 79.5 ± 13.7 cm measured in September. The water table in the peat plateau followed the deepening of the seasonally thawed layer.

Soil temperatures followed the seasonal climate but were dampened and had temporal lags in deeper peat layers (Fig. S1a). The highest young bog and mature bog soil temperatures at 10 cm depth occurred in July at 14.3 and 14.1 ^∘C, respectively. At 100 cm depth the maximum temperatures occurred in August and September at 8.6 and 6.9 ^∘C, respectively, for the young and mature bogs. Soil temperatures at 250 cm were still rising at the end of September, peaking at 4.1 and 3.2 ^∘C in the young bog and mature bog, respectively. The young bog was consistently warmer than the mature bog throughout the study by an average of 0.9 ± 0.9, 1.8 ± 1.0, and 0.5 ± 0.4 ^∘C at 10, 100, and 250 cm depths, respectively.

Across all depths and sampling occasions, average pH was higher (ANOVA: $F_{(\mathrm{1},\mathrm{77})}=\mathrm{35.2}$, P<0.001) in the young bog than in the mature bog at 4.1 ± 0.2 and 3.9 ± 0.2, respectively. In contrast, DOC at 69.2 ± 18.4 and 53.8 ± 5.4 mg C L⁻¹ (ANOVA: $F_{(\mathrm{1},\mathrm{82})}=\mathrm{38.7}$, P<0.001) and total dissolved nitrogen at 1.5 ± 1.4 and 0.9 ± 0.1 mg L⁻¹ (ANOVA: $F_{(\mathrm{1},\mathrm{82})}=\mathrm{12.8}$, P<0.01) were higher in the mature bog than in the young bog, respectively. Average SUVA values were higher (ANOVA: $F_{(\mathrm{1},\mathrm{82})}=\mathrm{103.5}$, P<0.001) in the young bog (3.2 ± 0.4 L mg C⁻¹ m⁻¹) compared to the mature bog (2.6 ± 0.4 L mg C⁻¹ m⁻¹), indicating DOM with a greater aromatic content in the young bog. However, average spectral slope (S_250−465) values were also greater (ANOVA: $F_{(\mathrm{1},\mathrm{81})}=\mathrm{6.9}$, P<0.05) in the young bog (−0.016 ± 0.002 nm⁻¹) compared to the mature bog (−0.017 ± 0.003 nm⁻¹), indicating lower molecular weight and decreasing aromaticity. Average phenolics (0.6 ± 0.2 and 0.6 ± 0.2 mg L⁻¹) and phosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-}$: 9.0 ± 14.3 and 6.7 ± 3.0 µg L⁻¹) were similar between the young bog and mature bog, respectively, across all depths and sampling occasions. Full details of DOM chemistry results can be found in Heffernan et al. (2021). Of note is the fact that the pore water chemistry was compared across all depths in this study, in contrast to Heffernan et al. (2021) in which pore water found above and below the transition indicating permafrost thaw was compared.

3.2 Concentrations and isotopic signatures of dissolved gases

Dissolved CH₄ increased with depth below the water table in both the young and mature bogs (Fig. 2a). Dissolved CH₄ concentrations in the young bog increased with depth, from 19 µmol L⁻¹ at 5 cm depth to a peak of 5400 µmol L⁻¹ at 195 cm. Dissolved CH₄ concentrations in the mature bog remained low above the water table (<6 µmol L⁻¹ below 25 cm) but then increased to 4100 ± 1700 µmol L⁻¹ between 115 and 250 cm depth and peaked at 6800 µmol L⁻¹. Dissolved CO₂ concentrations followed a very similar pattern to CH₄, increasing with depth in both the young and mature bogs (Fig. 2b). Again, the mature bog had overall higher concentrations, with mean average values ranging from 340 to 1295 µmol L⁻¹ and peaking at 1500 µmol L⁻¹ at 85 cm. In contrast, the young bog had average values ranging from 113 to 960 µmol L⁻¹ and peaked at 1200 µmol L⁻¹ at 95 cm (Fig. 2b).

Figure 2 Average seasonal (May–September) depth profiles in the young (green, black circles) and mature (yellow, black circles) bogs of (a) dissolved CH₄ concentration ( µmol L⁻¹), (b) dissolved CO₂ concentration (µmol L⁻¹), (c) δ¹³C-CH₄ ( ‰), (d) δ¹³C-CO₂ ( ‰), and (e) apparent fractionation factor (α_c) between dissolved CH₄ and CO₂. (f) Cross-plot of corresponding δ¹³C-CH₄ and δ¹³C-CO₂ values (‰) in the young bog and mature bog from raw data used in panels (c) and (d). Diagonal lines represent different α_c, where α_c 1.040–1.065 represents acetoclastic methanogenesis, and α_c 1.055–1.09 represents hydrogenotrophic methanogenesis (Whiticar, 1999). (a–e) Dashed and dotted horizontal lines represent water table depth in the young bog (YB) and mature bog (MB), respectively. Arrows in panel (a) represent depth of thaw transition in both the young bog (29 cm) and mature bog (71 cm), i.e., the transition from deep peat (accumulated prior to thawing) and shallow peat (accumulated post-thawing).

The young bog and mature bog had distinct profiles of δ¹³C values for both CH₄ and CO₂ (Fig. 2c, d). The young bog had no apparent trend with depth for both δ¹³C-CH₄ (ANOVA: $F_{(\mathrm{14},\mathrm{45})}=\mathrm{1.75}$, P=0.08) and δ¹³C-CO₂ (ANOVA: $F_{(\mathrm{14},\mathrm{46})}=\mathrm{1.79}$, P=0.07), averaging −62.4 ± 7.0 ‰ and −6.8 ± 1.6 ‰, respectively (Fig. 2c, d). In the mature bog we observed significant depth trends for both δ¹³C-CH₄ (ANOVA: $F_{(\mathrm{14},\mathrm{43})}=\mathrm{3.19}$, P<0.01) and δ¹³C-CO₂ (ANOVA: $F_{(\mathrm{14},\mathrm{49})}=\mathrm{6.22}$, P<0.001). These significant depth trends are due to isotopically heavy δ¹³C-CH₄ and light δ¹³C-CO₂ above the water table, which suggests an influence from CH₄ oxidation. When comparing δ¹³C depth profiles between the thermokarst bogs we focused on those values taken from under the water table to avoid the effect of CH₄ oxidation observed above the water table in the mature bog. Under the water table, δ¹³C-CH₄ values in the mature bog were significantly lighter (ANOVA: $F_{(\mathrm{1},\mathrm{64})}=\mathrm{18.72}$, P<0.001) compared to the young bog at an average of −68.7 ± 5.0 ‰ and −62.4 ± 7.0 ‰, respectively. Conversely, the mature bog had isotopically heavier δ¹³C-CO₂ than the young bog below the water table (ANOVA: $F_{(\mathrm{1},\mathrm{71})}=\mathrm{13.86}$, P<0.001).

The apparent fractionation factor (α_C) is a robust parameter to characterize the relative contribution of CH₄ production pathways, with values of 1.040–1.060 indicating acetoclastic methanogenesis and 1.060–1.090 indicating hydrogenotrophic methanogenesis (Chanton et al., 2005). Similar to the gas δ¹³C depth profiles, we found no clear trend with depth for α_C values in the young bog (ANOVA: $F_{(\mathrm{14},\mathrm{44})}=\mathrm{0.87}$, P=0.59) with an average of 1.058 ± 0.012 and range of 1.018–1.079 (Fig. 2e). In the mature bog, we found a clear depth trend in α_C values (ANOVA: $F_{(\mathrm{14},\mathrm{43})}=\mathrm{5.71}$, P<0.001). Similar to the δ¹³C depth profiles in the mature bog, this significant depth trend in α_C is due to the influence of CH₄ oxidation above the water table, with the lowest α_C values being those from samples collected above the water table at 5, 15, and 25 cm. The average α_C beneath the water table in the mature bog was 1.064 ± 0.017 and ranged from 1.015 to 1.094. When comparing α_C values from beneath the water table between the young and mature bogs, we found that α_C values were significantly lower in the young bog (ANOVA: $F_{(\mathrm{1},\mathrm{63})}=\mathrm{30.8}$, P<0.001).

In the isotopic ratio cross-plot of δ¹³C-CH₄ and δ¹³C-CO₂ (Fig. 2f), most of the young bog had α_C values between 1.055 and 1.065 (29 in total), with a greater number of samples (21) between α_C=1.040–1.055 compared to the mature bog (15). In contrast, a greater proportion of the mature bog samples had α_C>1.065 (42 in the young bog and 52 in the mature bog). There was no clear depth trend in the α_C values, and no samples in this study had α_C>1.090. Several samples (13) from the young bog and mature bog had α_C values of <1.040 likely due to CH₄ oxidation (Knorr et al., 2009).

3.3 Magnitude and isotopic signature of land–atmosphere gas fluxes

The young bog had almost 3 times greater average CH₄ fluxes than the mature bog during the May–September study period, at 82.3 ± 21.9 mg CH₄ m⁻² d⁻¹ and 30.8 ± 10.6 mg CH₄ m⁻² d⁻¹, respectively (Fig. 3a). Fluxes of CH₄ in the young bog were greatest between June and August, ranging from 80.6 ± 40.3 m CH₄ m⁻² d⁻¹ to 100.9 ± 63.1 mg CH₄ m⁻² d⁻¹. The lowest young bog CH₄ fluxes were observed in September at 55.0 ± 17.7 mg CH₄ m⁻² d⁻¹ (Fig. S3a). Mature bog CH₄ fluxes were greatest in September (55.8 ± 21.1 mg CH₄ m⁻² d⁻1) and lowest in May (5.6 ± 2.7 mg CH₄ m⁻² d⁻¹). Ecosystem respiration (CO₂ emissions measured with dark chambers) was significantly lower in the young bog than mature bog, with study period averages of 0.6 ± 0.3 and 1.9 ± 0.3 gĊO₂ m⁻² d⁻¹, respectively (Fig. S3). Maximum ecosystem respiration in the young bog occurred in August (1.6 g CO₂ m⁻² d⁻¹) and was much lower during the other 4 months (monthly averages of 0.2 to 0.4 g CO₂ m⁻² d⁻¹). Ecosystem respiration rates in the mature bog were elevated from June to August (monthly averages between 2.1 and 2.6 g CO₂ m⁻² d⁻¹) and decreased in September (0.8 g CO₂ m⁻² d⁻¹). The proportion of total C emissions (sum of CH₄ and CO₂ emissions) released as CH₄ were an order of magnitude greater in the young bog than mature bog stage, at 18 % and 2 %, respectively. This was a result of both higher CH₄ emissions and lower ecosystem respiration (Fig. S3) in the young bog. The δ¹³C-CH₄ signature of CH₄ emissions (intercept values from Keeling plots) in the young bog were significantly greater than those observed in the mature bog (Fig. 3c; ANOVA: $F_{(\mathrm{1},\mathrm{4})}=\mathrm{20.67}$, P<0.05). The average δ¹³C-CH₄ signature of CH₄ emissions in the young bog (n=4) was −66.5 ± 1.4 ‰ (95 % confidence interval) and 78.5 ± 5.6 ‰ (95 % confidence interval; Fig. 3c) in the mature bog emissions (n=4).

Figure 3 Magnitude and isotopic signature of greenhouse gas fluxes from the young bog (green) and mature bog (yellow) shown as boxplots. Boxes represent the interquartile range (25 %–75 %), with median shown as black horizontal line. Whiskers extend to 1.5 times the interquartile range (distance between first and third quartile) in each direction, with outlier data plotted individually as black dots. (a) The magnitude of net land–atmosphere CH₄ emissions as measured by soil chambers. (b) The ratio between CH₄ emissions and the sum of CO₂ emissions (ecosystem respiration) and CH₄, both standardized to per gram C. (c) Intercept values of Keeling plots indicating the δ¹³C-CH₄ signature of CH₄ emissions. Isotopically heavier (i.e., less negative) δ¹³C-CH₄ is produced via acetoclastic methanogenesis, whereas isotopically lighter (i.e., more negative) δ¹³C-CH₄ is produced via hydrogenotrophic methanogenesis, The CH₄ and CO₂ land–atmosphere fluxes shown in (a) and (b) were measured once a month from May to September 2018. The δ¹³C-CH₄ of CH₄ emitted to the atmosphere was measured in September and October 2016 (see “Methods” for details and Fig. S4 for Keeling plots).

3.4 Microbial community structure along the permafrost peatland thaw gradient

We used NMDS ordinations to assess differences in microbial community structure between solid peat and pore water samples, between sampling depths, and between the plateau, young bog, and mature bog. The only exception was the plateau, where only peat samples were collected (i.e., no pore water samples). Microbial community structure in peat was determined to be significantly different from porewater microbial communities (PERMANOVA: R²=0.13, P<0.05, Fig. 4). The differences observed in the microbial community structure between peat and pore water samples could be a function of the different extraction methods used to extract DNA (Carrigg et al., 2007). Among the pore water samples, distinct microbial communities were found to be associated with the young bog and mature bog. Similarly, microbial community structure in peat was found to be significantly distinct between the three successional stages (plateau peat, young bog, and mature bog; Fig. 4; PERMANOVA: R²=0.18, P<0.05). There is also a common trend in vertical community structuring for all sample matrices according to depth. Changes in overall microbial community composition in both peat and pore water, across a vertical profile (to a maximum depth of 160 cm), illustrate a confluence in microbial community structure with depth in both the young and mature bogs (Fig. 4). In other words, community structure was most dissimilar at depths closer to the surface (Figs. 4, S2; PERMANOVA: R²=0.16, P<0.05). This trend was particularly evident in the porewater samples (Fig. 4). In the peat samples, though microbial communities did not fully converge, deeper young bog peat (i.e., 90–160 cm) communities did become more similar to communities found in the mature bog at intermediate depths (i.e., 30–70 cm) based on the nearness of sample points in NMDS ordinations (Fig. 4). We also observed that the mature bog near-surface peat samples were located closer to the plateau peat in NMDS ordinations (Fig. 4; PERMANOVA: R²=0.4, P=0.1). It was not possible to assess the presence of this cyclic succession (from young bog to mature bog to plateau) in the pore water samples since we did not characterize the microbial community in the plateau pore water. Finally, we also assessed the effect of seasonality on microbial community structure and found no effect with regards to sampling month (PERMANOVA: R²=0.02, P=0.090).

Figure 4 Microbial community distribution according to stage of peat/pore water. NMDS ordinations of amplicon sequencing variant (ASV) data demonstrate significant community dissimilarities according to thaw stage for both pore water (shown by the triangles) and peat (shown by the circles) samples, encircled by 95 % confidence intervals. Colour gradient and lines demonstrate the shift in microbial community structure along vertical depth profiles in which lighter shades indicate samples closer to the surface.

The total archaeal community comprised 6 % of the entire microbial dataset. Methanogen-related orders comprised 54 % of this archaeal dataset and demonstrated marked differences in the relative abundance of acetoclastic-related methanogens according to thaw stage and depth in both peat and pore water samples (Figs. 5, S2). In the young and mature bog peat samples, hydrogenotrophic-related methanogens were ubiquitously present throughout both depth profiles (Fig. 5a). In comparison, acetoclastic-related methanogens exhibited a relatively restricted presence and were only present at specific depths (Fig. 5a). These communities were most abundant (>25 % of the total methanogenic community) near the surface in the young bog, just above and below the thaw transition zone (Fig. 5a). In the pore water, hydrogenotrophic methanogens were also dominant throughout depths in both stages of thaw (Fig. 5b). However, in contrast to peat samples, acetoclastic methanogens were virtually absent in the pore water, although minimally present (i.e., ≤ 10 % relative abundance) at depths between 35 and 155 cm, all found below the thaw transition zone (Fig. 5b).

Figure 5 Relative abundance of archaeal orders according to putative methanogenic capability along a depth profile for peat and pore water samples. Samples are arranged according to depth (y axis), with the relative abundance of methanogenic archaea resolved shown on the x axis. Note that the y axis does not uniformly progress in 10 cm increments. Progression of thaw is shown from plateau peat to young bog to mature bog at the top of the figures, with position of water table shown in blue for each panel. Red lines demonstrate thaw transition zone for the young bog and mature bog. (a) Stacked bar plot of methanogenic Archaea for all peat samples. Samples demonstrate significant differences in putative methanogen composition between all stages (Kruskal–Wallis test and Wilcox rank sum test, with Benjamini–Hochberg-corrected P values, P<0.05). (b) Stacked bar plot of methanogenic Archaea for all pore water samples. Samples do not demonstrate significant differences in putative methanogen composition between stages (Kruskal–Wallis test, with Benjamini–Hochberg-corrected P values, P=0.965).

Using a redundancy analysis (RDA, Fig. 6) we found that 27.6 % of variation in the methanogenic community was explained by two variables: thaw stage (ANOVA: P<0.05) and depth from the water table (ANOVA: P<0.05). Although these were the only two parameters that were identified as significant variables impacting microbial community structure when using a backward stepping model, it should be noted that there may be more variation in the community that our experimental design does not take into account as a result of unconstrained variation represented by plant–microbe and/or microbe–microbe interactions (Boon et al., 2014). Nonetheless, the 27.6 % variation explained is in accordance with other studies conducted in permafrost-impacted regions using similar methods, for which the percentage of explained variation falls between 6 % (low) to 43 % (high) (Comte et al., 2015; Hough et al., 2020). Next, we used variance partitioning to assess the extent to which thaw stage and depth from the water table (i.e., the significant environmental variables identified by the RDA) explained the variation in only the methanogenic community structure (Fig. 6). Based on this analysis, thaw stage explained 18.4 % and distance to the water table explained 4.3 % of methanogenic community variation.

Figure 6 Redundancy analysis (RDA) exploring significant biotic and abiotic variables influencing the total methanogenic community (adjusted R²=27.6 %), as determined by a backward stepping RDA model in the peat and pore water samples. All parameters that were used in the model are described in Sect. 2.1 of the “Methods”. Grey dots in the panel demonstrate the distribution of all ASVs in the methanogenic dataset. Shaded ellipses represent the 95 % confidence intervals for microbial community structure according to peatland thaw stage (young bog versus mature bog). Only significant (ANOVA: P<0.05) variables are shown. Using variation partitioning, we found that peatland thaw stage significantly explains about 18.4 % of methanogenic community variation, whereas distance to water table explained 4.3 %. Both axes are significant (ANOVA: P<0.05).

4 Discussion

Our study shows that high CH₄ emissions from thermokarst bogs in the initial decades following permafrost thaw (young bog) are not only linked to environmental conditions (wetness, soil temperature, vegetation) but are also driven by relatively increased microbial CH₄ production through the energetically more favourable acetoclastic methanogenesis pathway. Evidence of putatively acetoclastic methanogens and CH₄ produced via the acetoclastic metabolic pathway was found in the young bog both near the surface and at depths below the thaw transition (i.e., in peat that accumulated prior to permafrost thaw). We are unable to determine whether these greater CH₄ emissions in the initial decades following thaw are due to the mineralization of labile organic matter released from previously frozen peat or are driven solely by fresh, labile DOM derived from surface vegetation leached throughout the peat profile. However, previous work in the discontinuous permafrost region in the Interior Plains of western Canada has found a limited contribution of previously frozen organic matter contributing to surface CH₄ emissions in thermokarst bogs (Cooper et al., 2017). Elevated CH₄ emissions then slow over the following centuries with succession into a mature thermokarst bog stage in which CH₄ production is almost exclusively through the hydrogenotrophic pathway.

4.1 Shift in microbial community assemblages along a permafrost thaw gradient

Microbial communities varied along the permafrost thaw gradient: among different thaw stages (permafrost peat plateau, young bog, and mature bog), with peat depth (surface down to 160 cm), and between different sample types (solid peat and pore water). We found clear differences in microbial communities between the young bog and mature bog despite similar peat stratigraphy up to the surficial vegetation (Heffernan et al., 2020), where dominant Sphagnum species varied. The greater height of the peat surface above the water table and drier conditions in the mature bog, due to the slow accumulation of new peat over centuries, lead to a shift in vegetation composition from hydrophilic Sphagnum and graminoids towards more drought-resistant Sphagnum spp. and ericaceous shrubs. This shift in water table position and vegetation community, along with a decrease in temperatures (Fig. S1a) due to the thermal insulating properties of Sphagnum peat (Kujala et al., 2008), appears to have caused the observed differences in microbial communities between the young and mature bogs, even at depths >1 m. Microbial communities were most dissimilar between the peat plateau and young bog. This was unsurprising given the abrupt shift from the elevated, frozen, and relatively dry peat plateau forest to the young bog where the surface was saturated, was dominated by hydrophilic vegetation, and had warmer temperatures. We further noted that the microbial community of the mature bog was more similar with the peat plateau than with the young bog. Paleo-records in the region (Heffernan et al., 2020; Pelletier et al., 2017; Zoltai, 1993) show that many peatlands have undergone cyclical permafrost developments as thermal insulating properties of Sphagnum peat in mature bogs leads to the re-aggradation of permafrost peat plateaus. Our study suggests that the peat plateau microbial community is influenced by the preceding mature bog microbial community as permafrost aggrades.

The most dissimilar microbial community composition was observed between samples near the surface and those at depth (i.e., down to 160 cm), as has also been observed in other permafrost ecosystems (Frey et al., 2016; Monteux et al., 2018). Shifts in microbial community composition along the thaw gradient were most evident nearer the surface, whereas communities found at depth were similar between the young bog and mature bog (Fig. 4). At the surface, microbial community structure is influenced by the successional vegetation community (Hodgkins et al., 2014), and the role that vegetation, particularly graminoids which are found in the young bog, has on microbial community structure has been well documented in northern peatlands (Robroek et al., 2015, 2021; Bragazza et al., 2015). Moderately acidic, saturated peatlands with hydrophilic vegetation, similar to the young bog, have been shown to harbour acid-tolerant fermenting bacteria that produce substrates for methanogenesis and are trophically linked with methanogens (Wüst et al., 2009). Thus, the interaction between water table position, pH, and vegetation community influences the substrates available to the microbial community, which in turn impacts the surface community's structure (Kotiaho et al., 2013). In contrast, communities at depth are known to be influenced by peat properties, such as peat chemistry and degree of decomposition, and the paleoenvironment under which they originally colonized (Lee et al., 2012; Holm et al., 2020). In the young and mature bogs, both peat properties (humification indices including FTIR (Fourier-transform infrared spectroscopy) $\mathrm{1630}/\mathrm{1090}$ cm⁻³ and C : N ratios) and the paleoenvironment at depth are similar (Heffernan et al., 2020), which may explain the observed convergence of microbial community structure. Nonetheless, although there are some similarities at depth between both young and mature bogs, microbial communities inhabiting either are still distinct (Fig. 4). This is emphasized by the differing abundance of Archaea that participate in hydrogenotrophic or acetoclastic methanogenesis (Fig. 5) in both stages down the peat profile.

As has been shown previously in other thermokarst peatlands (McCalley et al., 2014), the young and mature bog stages were dominated by hydrogenotrophic methanogens. However, putatively acetoclastic methanogens were relatively more abundant in the young bog (Fig. 5), particularly at or below the transition in peat that accumulated prior to permafrost thaw. Thaw stage and distance from the water table were found to influence the methanogenic community composition (Fig. 6), with distance from the water table dictating where anoxic conditions persist (Blodau et al., 2004) and thus where methanogenic colonization can occur. The influence of vegetation communities associated with different thermokarst peatland stages on methanogenic community composition has previously been attributed to the role of plant-derived DOM serving as the substrate for CH₄ production (Liebner et al., 2015; McCalley et al., 2014). The presence of hydrophilic vegetation, particularly graminoids, in the saturated young bog provides the precursors for fermentation, yielding acetate (Liebner et al., 2015; Ström et al., 2003, 2012, 2015) and serving as the substrate for acetoclastic CH₄ production. The downward transport from the surface of plant-derived DOM in the young bog (Chanton et al., 2008) likely provides sufficient acetate for the establishment of acetoclastic methanogens at depth in this environment.

4.2 Production and emissions of CH₄ along a peatland thaw gradient

Isotopic signatures (δ¹³C) of dissolved CO₂ and CH₄, α_C values in porewater, and the δ¹³C signature of CH₄ emitted to the atmosphere provided further evidence of relatively elevated acetoclastic methanogenesis in the young bog stage. The general increase in δ¹³C-CO₂ with depth observed at both sites (Fig. 2d) indicates accumulation of isotopically heavier δ¹³C-CO₂ which is likely explained by the preferential use of isotopically lighter δ¹³C-CO₂ during hydrogenotrophic methanogenesis (Hornibrook et al., 2000). As a result, CH₄ tends to become lighter with depth, and this was particularly apparent in the mature bog (Fig. 2c). This leads to the average α_C values of 1.064 (δ¹³C-CH₄; −68.7 ‰) in the mature bog, which were significantly higher than the 1.058 (δ¹³C-CH₄; −62.4 ‰) observed in the young. Together, the δ¹³C-CH₄ and δ¹³C-CO₂ data and the resulting α_C depth profiles suggest that the majority of CH₄ is produced via the hydrogenotrophic methanogenic pathway, which supports the findings of the microbial community analysis (Fig. 5). Our isotope data also suggest that a greater proportion of CH₄ is produced via acetoclastic methanogenesis throughout the profile in the young bog compared to the mature bog (Fig. 2c–f). This is evident from lower average α_C values found in the young bog compared to the mature bog, as well as the greater number of these young bog α_C values falling between 1.040 and 1.065, which represents acetoclastic methanogenesis (Whiticar, 1999). These findings again agree with the relatively greater abundance of acetoclastic methanogens observed at that site (Fig. 5).

In this study we found that average CH₄ emissions in the initial decades following thaw, in the young bog stage, were 2.5–3 times greater than emissions measured in the mature bog stage which had thawed ∼ 200 years ago (Fig. 3a). Furthermore, the proportion of CH₄ to overall C emissions (Fig. 3b) was considerably greater in the young bog than in the mature bog. In the mature bog the lower water table position leads to both increased CO₂ emissions and decreased CH₄ emissions, resulting in a reduced fraction of C emissions as CH₄. Previous studies have shown similarly increased CH₄ emissions in the initial decades following thaw (Johnston et al., 2014; Wickland et al., 2006). While our pore water chemistry data are inconclusive with regards to organic carbon characteristics, other work in thermokarst bogs in the Interior Plains of western Canada has shown that the organic matter derived from the young bog vegetation community is highly labile (Burd et al., 2020). Previous work at our study site has shown that the vegetation community in the young bog is associated with greater potential enzymatic degradation of organic matter (Heffernan et al., 2021). Hydrolysis of plant-derived organic matter by extracellular enzymes leads to the formation of monomers (Kotsyurbenko, 2005). These monomers can be further degraded to form acetate and other precursors for methanogenesis when present with anaerobic fermenting bacteria (Hamberger et al., 2008) and near the surface and vegetation inputs (Hädrich et al., 2012). Our study shows that these higher CH₄ emissions are likely linked to increased wetness, temperatures, and a vegetation community associated with more labile organic matter which favour a greater proportion of CH₄ produced via acetoclastic methanogenesis, as shown by our δ¹³C-CH₄ values, α_c depth profiles, and microbial community composition analyses.

Many factors, including environmental conditions and microbial community structure likely contribute to the differences in net CH₄ emissions from the young and mature bogs (Fig. 3a). Methane oxidation has been shown to be an important regulator of post-thaw CH₄ emissions (Perryman et al., 2020) and to result in isotopically heavier (i.e., less negative) δ¹³C-CH₄ and lighter (i.e., more negative) δ¹³C-CO₂ (Whiticar, 1999). Our data suggest the role of CH₄ oxidation was different between sites. Methane oxidation was apparent in the δ¹³C-CH₄ and δ¹³C-CO₂ signatures above the water table in the mature bog, but no CH₄ oxidation is evident in the young bog (Fig. 2c, d). The difference in gas flux δ¹³C signatures (Fig. 3c) also suggests a greater prevalence of CH₄ oxidation in the mature bog. However, increased oxidation above the water table in the mature bog is likely not fully responsible for the observed differences in CH₄ surface emissions and depth profiles between the young and mature bogs. Lower soil temperatures, a vegetation community associated with reduced substrate availability, the dominance of hydrogenotrophic methanogenesis throughout the peat profile, and a deeper water table position all contribute to the lower CH₄ production and higher CH₄ oxidation observed in the mature bog. Nonetheless, using this interdisciplinary approach, we are unable to determine the relative contribution of acetoclastic methanogenesis at each depth to the overall emissions at the surface.

Our results, and those of others (Euskirchen et al., 2014; Johnston et al., 2014), have shown that CH₄ emissions exhibit seasonal variation (Fig. S3a, c). However, in contrast to some previous findings (Ebrahimi and Or, 2017), we did not observe a corresponding seasonal response in the microbial community composition (Fig. S2). This may be a sampling design effect since our study spanned only 2 months (June and September), compounded by the fact that we did not have replicate samples to test the robustness of this finding. However, other studies have also shown that soil microbial community growth is not impacted by seasonal variations in temperature (Simon et al., 2020) and that microbial communities require a longer timescale (years–decades–centuries) to respond to temperature following thaw (Feng et al., 2020). Our results corroborate these observations, suggesting a long-term response in the microbial community composition to the ecological shifts associated with autogenic peatland succession following permafrost thaw. Autogenic peatland succession following thaw occurs on the decade to century timescale, shifting from recently thawed to mature thermokarst bogs (Camill, 1999). Both recently thawed (young) and mature thermokarst bogs have distinct hydrological regimes, vegetation communities, and peat chemistry. Following thaw, associated changes in vegetation and litter input alter microbial community composition and activity (Adamczyk et al., 2020; Kirkwood et al., 2021). Such changes in microbial community structure thus impact CH₄ emissions from thermokarst peatlands. Under predicted climatic warming scenarios differences in microbial community composition have been shown to be increasingly driven by seasonally independent variables such as substrate quality and the legacy effects of soil temperatures (Luláková et al., 2019). This study suggests that the environmental conditions required for increased methanogenic activity at depth is limited to the initial decades following thaw, after which the microbial community structure changes in response to the lowering of the water table, lower soil temperatures, and shifts in the vegetation community.

5 Conclusions

This study demonstrates that higher CH₄ emissions in thermokarst bogs in the initial decades following thaw are driven by shifts in vegetation communities that produce organic matter inputs of varying lability (Burd et al., 2020) and a prevalence of anoxic conditions, which was associated with an increase in acetoclastic methanogenesis in our site. The influence of this pathway was apparent at depth throughout the peat profile. With succession following thaw towards a mature thermokarst bog, a shift in water table position and vegetation composition seems to reduce the role of the acetoclastic methanogenesis pathway. Previous work at this site (Heffernan et al., 2021) and other thermokarst peatlands in the discontinuous permafrost zone of boreal western Canada (Burd et al., 2020) have indicated that the vegetation community found in the initial decades following permafrost thaw is associated with increased potential enzymatic degradation and biodegradability of organic matter compared to that found in the mature bog. Average growing season CH₄ emissions were 2.5–3 times greater in the recently thawed young bog. Overall, C emissions in the young bog contained proportionally more CH₄ than those from the mature bog due to greater CH₄ production and also reduced CO₂ emissions. These greater CH₄ emissions in the young bog are driven by a higher contribution to surface emissions from CH₄ produced throughout the peat profile by acetoclastic methanogens. The response of the microbial community to permafrost thaw is tied to the shifting environmental conditions associated with peatland autogenic succession. Warmer and wetter conditions in the initial decades following thaw, in conjunction with a vegetation community associated with greater availability of labile plant leachates (Bragazza et al., 2015), provide favourable conditions for acetoclastic methanogens throughout the peat profile. Given the projected increases in thermokarst peatland formation (Olefeldt et al., 2016), our study suggests that we can expect a pulse of CH₄ emissions from current regions of the discontinuous permafrost zone. This pulse will be driven, in part, by increased acetoclastic methanogenesis from labile substrates in recently thawed thermokarst peatlands. However, this rapid increase in CH₄ emissions will only remain at the decadal to century scale as autogenic peatland succession results in relatively drier mature thermokarst bogs, where lower temperatures and less labile substrate availability lead to a dominance of hydrogenotrophic methanogenesis.