Stordalen Mire (68.35°N, 19.05°E) is a peat plateau in subarctic Sweden underlain by discontinuous permafrost, which is thawing as the region warms. The Mire has three dominant vegetation community types (Johansson et al., 2006): (1) Palsa (intact permafrost) with woody herbaceous vegetation; (2) Ombrotrophic peatland or bog (intermediate thaw) features with Sphagnum spp., sedges, and shrubs as the dominant vegetations; and (3) Minerotrophic peatland or fen (fully thawed) dominated by sedges, mainly Eriophorum spp. (Bäckstrand et al., 2010). Generally, bogs are fed only by precipitation, have acidic (pH ~ 4) surface waters and Sphagnum-dominated vegetation. These Sphagnum mosses, are known to inhibit alternate plant growth by acidifying water and creating an environment that is harsh to other plant species (van Breemen, 1995). Fens, on the other hand, are characterized by more mineral-rich water because they are fed by surface water and groundwater, in addition to rainfall. They are usually dominated by grasses and sedges and maintain slightly acidic to alkaline pH. Due to continued warming in the Arctic, the areal extent of intact palsa (Malmer et al., 2005) and Sphagnum community (Norby et al., 2019) across the Stordalen Mire has declined significantly since 1970 (Malmer et al., 2005) while at the same time, wetter fen habitats have expanded (Kokfelt et al., 2009; Bäckstrand et al., 2010). The ground's thaw state and consequent relationship to the water table determine the plant community, which in turn determines organic matter composition. In areas where the permafrost is degraded, land subsidence leads to varying degrees of water inundation. In this study, we focused on these inundated sites, and specifically a fen and a bog site with an extensive set of historical gas flux, microbial, and geochemical measurements (previously described by Hodgkins et al., 2014, 2016; McCalley et al., 2014; Woodcroft et al., 2018). In June 2012, peat was collected from both sites with an 11-cm-diameter homemade circular push corer. Both sites were cored in triplicate to capture spatial heterogeneity. Cores were divided into three sections, and the sections were placed into plastic bags and stored frozen until analysis. Eighteen soil samples were collected in total, including nine samples (3 cores × 3 depths per core) collected from a Sphagnum-derived ombrotrophic bog and nine samples collected from a sedge-dominated fen (Table 1). Porewater for dissolved organic carbon (DOC) analysis was gathered from the same depths by syringe suction through a 1-m long, 0.5-cm-diameter stainless steel tube with holes drilled along the bottom 3 cm, then filtered through 0.7-μm Whatman GF/F glass microfiber filters into 120-mL brown borosilicate bottles. The pH was measured on-site with an Oakton Waterproof pHTestr 10. The bottles were frozen within 8 hours of collection and stored frozen until analysis. DOC concentrations were measured by high-temperature catalytic oxidation on a Shimadzu Total Organic Carbon analyzer with a non-dispersive infrared detector. Each sample was analyzed in triplicate measurements, which constantly had a coefficient of variance of <2%. The concentration was reported as millimolar (mM). Sulfate and nitrate concentrations were measured by ion chromatography on a Dionex ICS-1100 with a 4-mm IonPac AS22 column, with a 4.5 mM carbonate/1.4 mM bicarbonate eluent at a flow rate of 1.2 mL/min. The concentration was reported as micromolar (μM). In this study, we will focus on the overall SOM molecular composition differences observed in the bog vs. the fen and not the depth effect.

Porewater samples for dissolved gas analysis were collected similarly to the DOC samples, with porewater drawn up by syringe suction through perforated stainless-steel tubes that were inserted into the peat at the desired depth. Porewater was immediately filtered in the field using 2.7-μm Whatman GF/D glass-fiber filters, then stored in pre-evacuated glass vials sealed with butyl stoppers. Samples were then processed according to the method described by (Hodgkins et al., 2015). Phosphoric acid (1 mL of 20%) was added to each sample, to convert all dissolved inorganic carbon to CO₂ (thus removing differences in HCO3-/CO₂ ratios due to differing pH) and to preserve for shipment to Florida State University. Helium was then added to the sample headspaces to bring them to a constant pressure of 1 atm. Samples were analyzed for CH₄ and CO₂ concentrations and stable isotopic composition (δ¹³C) on a Thermo Finnigan Delta-V Isotope Ratio Mass Spectrometer. As described in (Corbett et al., 2013, 2015) and (Hodgkins et al., 2014), headspace gas concentrations were then converted to dissolved concentrations (mM) in the original unacidified porewater based on previously-determined gas extraction efficiencies (defined as the proportion of formerly-dissolved gas in the headspace after acidification and equilibration between the water and headspace) (Supplementary Table 1). Each sample was analyzed twice and the average results for each sample were recorded. Analytical precision of replicate δ¹³C measurements was 0.2%0.

A 12 Tesla Bruker FT-ICR mass spectrometer located at PNNL was used to collect high resolution mass spectra of the WEOM. A Suwannee River Fulvic Acid standard (SRFA), obtained from the International Humic Substance Society (IHCC), was used as a control to monitor potential carry over from one sample to another and ensure instrument stability in day to day operation. The instrument was flushed between samples using a mixture of Milli-Q water and HPLC grade methanol. The ion accumulation time (IAT) was varied between 0.03 and 0.05 s to account for variations in C concentrations in different samples. For each sample, 144 individual scans were averaged, and an organic matter homologous series separated by 14 Da was used as an internal calibration. The mass measurement accuracy across 100–1,200 m/z range was <1 ppm (accounted for singly charged ions) and the mass resolution at 341 m/z was ~240 K and the transient was 0.8 s. BrukerDaltonik version 4.2 data analysis software was used to convert raw spectra (obtained from each sample) to a list of m/z values applying FT-MS peak picker module with a signal-to-noise ratio (S/N) threshold of 7 and absolute intensity threshold of 100. Formularity (Tolić et al., 2017) software was used to assign chemical formulae based on the criteria described in Tfaily et al. (2018).

The molecular formulae in each sample were evaluated on van Krevelen diagrams (Van Krevelen, 1950), based on their H:C ratios (y-axis) and O:C ratios (x-axis) (Kim et al., 2003), assigning them to the major biochemical classes (e.g., lipid, protein, lignin, carbohydrate etc.). The H:C and O:C ranges for biochemical classification are provided in Table 2.

Oxidation state is related to the bioenergetic potential of a specific compound based on its atomic makeup. Oxidation state indicates how a compound may be transformed in biochemical processes without necessarily knowing its structure (LaRowe and van Cappellen, 2011). For each individual compound, the average nominal oxidation state of carbon (NOSC), which represents the number of electrons transferred in C oxidation half reactions, was calculated based on Equation 1 (where C, H, N, O, P, and S represent the numbers of each element, and number coefficients represent typical oxidation states). More negative values for NOSC indicates the presence of more reduced carbon atoms.

ΔG°C-ox (kJ/mol C) indicates the thermodynamic favorability of the compound, with lower values of ΔG°C-ox for a compound indicating greater likelihood of degradation (LaRowe and van Cappellen, 2011). ΔG°C-ox was estimated from NOSC following Equation (2):

Positive values for ΔG°C-ox indicate that the oxidation of the C must be coupled to the reduction of a terminal electron acceptor, which in the TEA-depleted mire can include the soil organic matter itself (Keller and Takagi, 2013; Wilson et al., 2017). Organic matter degradation may become thermodynamically limited depending on the availability of terminal electron acceptors (Boye et al., 2017). Under oxic conditions the availability of high energy yielding O₂ as a terminal electron acceptor couples the decomposition of organic matter of all naturally occurring NOSC values (Keiluweit et al., 2016). However, under the inundated conditions of subsurface peat deposits, oxygen and other terminal electron acceptors become depleted thereby slowing or completely inhibiting the decomposition of organic matter with low NOSC values (Boye et al., 2017). In peat, oxygen is rapidly depleted to below detection limits (20 ppb) within 2–3 cm below the water surface (Lin et al., 2014). At Stordalen Mire, the average annual water table is typically near the peat surface in the fen (Hodgkins et al., 2014) and was 3.5–7 cm above the peat surface during our sampling. Thus, even the shallow fen peat is consistently under anoxic conditions and subject to depleted TEA availability. In the bog, the mean annual water table depth is around 17 cm (Hodgkins et al., 2014), but was only 2 cm below the peat surface during sampling for this project. Thus, the bog peat is subjected to fluctuating water (and subsequently oxygen) levels throughout the growing season. This could have the effect of enhancing decomposition rates during dry periods and re-oxidizing previously reduced TEAs to support decomposition following subsequent inundation. Relative energy yields associated with the terminal electron acceptors in anoxic conditions also provide insight on the sequences in which organic matter gets decomposed (Grandy and Neff, 2008). Lower ΔG°C-ox is an indication of compounds more thermodynamically favored to undergo metabolic and chemical reactions (Wilson and Tfaily, 2018) (i.e., available for microbial communities), whereas higher ΔG°C-ox reflects microbially degraded and transformed (i.e., recalcitrant) OM.

For each elemental formula, double-bond equivalent (DBE) minus oxygen (DBE-O) was calculated using Equation (3):

While NOSC and ΔG°C-ox provide information on overall quality of the OM, DBE-O provides more structural information, mostly of carbon-to-carbon double bonds (i.e., lack of hydrogen saturation) without considering the oxygen atoms that mostly exist as carboxyl groups in OM (Bae et al., 2011). Using this formulation, more oxygenated molecules have lower (negative) values for DBE-O while higher values for DBE-O indicate highly unsaturated compounds such as condensed aromatic structures (Herzsprung et al., 2014).

The rich database generated by high resolution mass spectrometry provides accurate mass information, enabling us to infer metabolite composition in complex matrices. We used network analysis to explore metabolic transformations in the bog and fen samples to further understand how two compounds are related to one another and how they might be involved in the same metabolic transformation (via a mass difference between two m/z values demonstrated as an edge in the network, and the initial and final compounds as nodes). Here, we used molecular transformation analysis to identify potential microbial transformation or decomposition pathways using FT-ICR-MS data. For a reaction to be considered and hence counted by this method, both the reactant and the product should be present, hence this approach can provide some information about SOM decomposition, uptake and accumulation (Burgess et al., 2017). We used Cytoscape (Shannon et al., 2003; Kohl et al., 2011) with the MetaNetter plug-in to generate molecular interaction networks. For each sample, the list of m/z values generated by the FT-ICR-MS analysis along with a list of 86 pre-defined m/z values and their subsequent molecular formula and IUPAC names [adopted from (Breitling et al., 2006) (Supplementary Table 2)] were imported into Cytoscape. The molecular interaction networks are consisting of nodes (m/z values) and edges (mass difference between two m/z values). Each m/z value inferred from the FT-ICR-MS analysis was defined as a node (described previously in Longnecker and Kujawinski, 2016) connected to other node (s) through the pre-defined m/z values defined as edges, representing a metabolic transformation. For example, gain or loss of CH₂ (Δm = 14.01565007) between two compounds in the molecular interaction networks is depicted by two nodes and an edge. Nodes represent the two compounds with mass difference of 14.01565007 Da and edge is the mass difference (named CH₂).

In total, 86 metabolic transformations were grouped into four different categories: amino acid-associated transformation (loss/gain of an amino acid and/or involvement of nucleoside or nucleotide), sugar -associated transformation (loss/gain of monomeric or polymeric carbohydrates), carbon transformation (reactions that can occur both biotically or abiotically depending on the environmental conditions), and other transformation that includes heteroatoms (such as loss/gain of a phosphate or an amine group). For this study, we have only focused on the first three categories since similar trends for the bog and fen samples (regardless of depth) were observed in the ‘other transformation’ category. The output of the Cytoscape is a list with the number of times (count) each transformation occurred between any two peaks in a particular sample as shown in Supplementary Figure S5. Normalized values (to the total number of transformations observed in each sample) and subsequently, normalized percentages are then calculated and used for comparison between samples and across sites.

Furthermore, a comprehensive set of topological parameters calculated from these undirected network plots (Assenov et al., 2008) based on transformation data generated by the Cytoscape with the MetaNetter plug-in was used to compare SOM metabolic transformation pathways between the bog and fen samples. For this study, two different parameters were chosen: clustering coefficient and average number of neighbors, all of which can be used to represent the level of population of the transformation network (Barabási and Albert, 1999; Ravasz et al., 2002). The clustering coefficient, C_n, in undirected networks is a ratio of the number of edges (N) between the neighbors of a node n, and the maximum number of edges (M) between the neighbors of n (Barabási and Oltvai, 2004). C_n is calculated as C_n = 2e_n/(k_n(k_n-1)), where k_n is the number (size) of neighbors of n and e_n is the number of connected pairs between all neighbors of n (Watts and Strogatz, 1998). The clustering coefficient of a node is always between 0 and 1(Barabási and Oltvai, 2004). A high clustering coefficient indicates a high number of neighbors for each node and hence represents low selectivity and high diversity since a compound could be produced by multiple reactions. The average number of neighbors therefore is an average of k_n of all nodes in the undirected network (Dong and Horvath, 2007). High number of neighbors indicates low degradation and high selectivity. Values for clustering coefficient and average number of neighbors for each sample are provided in Supplementary Table 3. Moreover, the network analysis parameters (clustering coefficient and average number of neighbors) were used to correlate the overall number and/or density of the metabolic reactions happening in each sample to measured biogeochemical properties in pore water such as pH, δ¹³C and DOC.

Principal component analysis of the relative abundance of different SOM compound classes (as inferred by FT-ICR-MS), biogeochemical properties of peat porewater (TN concentration, DOC, sulfate, nitrate, and ammonia concentration and pH) and peat porewater CO₂ and CH₄ concentrations revealed clear separation between the bog and fen sites along PC1 (which explained 41% of the variance), with higher PC1 scores in rich fens (Figure 1). SOM in the Sphagnum-dominated bog sites had a higher relative abundance of tannins, lignins, unsaturated, and condensed hydrocarbons (UnSatHC and CondHC, respectively), and lipids compared to the fen sites. In contrast, fen samples were separated from bog samples along PC1 which correlated positively with the abundance of peptide-, carbohydrates (cellulose-)-, and aminosugar-like (AS) compounds. PC2 explained a small proportion of the variance (16%). Higher concentrations (measured in mM) of TN and DOC were observed in the bog porewaters compared to the fen porewater. In contrast, sulfate, nitrate and ammonia concentrations were higher in fen porewater and correlated positively with peat porewater pH. Porewater CO₂ and CH₄ concentrations were higher in the fen sites. Although some depth trends were observed, these trends were generally weak, as they were not consistently observed over the entire depth range and none of the principal components showed significant correlations with depth (Supplementary Figure 1).

To better visualize the differences in the bog and fen SOM molecular composition, unique peaks present in each group (bog vs. fen, regardless of the depth) were reported separately on van Krevelen diagrams (Figures 2A,B). SOM in bog samples (Figure 2A) appeared to be more rich in lignin-like (0.29 < O/C < 0.65 and 0.7 < H/C < 1.5), tannin-like (0.65 < O/C < 1 and 0.8 < H/C < 1.5), and lipid-like (0 < O/C < 0.3 and 1.5 < H/C < 2.5) compounds and fen SOM (Figure 2B) was more rich in peptide-like (0.3 < O/C < 0.6 and 1.5 < H/C < 2.3) compounds and amino sugars (0.5 < O/C < 0.7 and 0.8 < H/C < 2.5). Out of 6,906 peaks detected in the bog samples in total, 24% were unique to bog, while only 8% of the peaks (out of 5,973 in total) detected in the fen were unique to fen, suggesting a higher number of unique compounds in bog (potentially not yet degraded) compared to fen samples (Venn diagram in panel B).

Figure 3A shows ΔG°C-ox calculated for 6906 assigned molecular formula in nine bog samples and 5973 assigned molecular formulas in nine fen samples. Lower ΔG°C-ox and high NOSC represent thermodynamically more labile organic compounds and vise-versa. For each sample, ΔG°C-ox values (Figure 3A; Table 3) and NOSC values (Table 3) were averaged across all formulas and unpaired t-test was performed on the averaged values obtained from nine bog samples and nine fen samples. Average ΔG°C-ox values in bog samples were lower than fen samples (Figure 3A; Table 3) and subsequently, average NOSC values were higher in bog samples compared to fen samples (Table 3). Interestingly, different ΔG°C-ox values were obtained for bog and fen samples within the same class of chemical compounds (Figure 3B). While amino sugars were more abundant in the fen sites compared to bog sites, they were characterized with higher average ΔG°C-ox compared to those in the bog sites. On the other hand, while bog sites had an abundance of lignin-like and tannin-like compounds compared to fen, however, they were of lower average ΔG°C-ox compared to those in the fen sites (more oxidized). The lower overall NOSC values of the bog SOM relative to the fen SOM could reflect the influence of oxygen at the bog site. As noted in section Chemical Compounds Classification and Data Processing, the fen is more consistently inundated (hence anoxic), while the bog habitat experiences alternating periods of oxic and anoxic conditions. These periods of exposure to oxygenation and the potential re-oxidation of organic and inorganic TEAs may facilitate decomposition at the bog site resulting in SOM with an overall lower NOSC value relative to the fen. However, the overall lower NOSC in the bog site may also reflect different starting substrates.

To further investigate the sulfur and nitrogen dynamics in the Mire, composition of S-containing SOM (regardless of depth) was compared in the bog (A) and fen (B) samples and unique peaks were reported on a van Krevelen diagram (Figures 4A,B). Sulfur containing compounds in bog fell in the lipid-like (0 < O/C < 0.3 and 1.5 < H/C < 2.5), peptide-like (0.3 < O/C < 0.6 and 1.5 < H/C < 2.3), and amino sugar-like (0.5 < O/C < 0.7 and 0.8 < H/C < 2.5) regions of the van Krevelen diagram. In contrast, S-containing compounds in fen fell in the condensed hydrocarbon-like (0 < O/C < 0.4 and 0.2 < H/C < 0.8) and lignin-like (0.29 < O/C < 0.65 and 0.7 < H/C < 1.5) regions of the van Krevelen diagram. Overall, a higher abundance of unique S-containing organic compounds was observed in the bog compared to the fen system (~ 2.6 to ~1.3% unique S containing compounds, respectively, depicted on Venn diagram in panel A) coincident with lower sulfate concentrations in bog porewater (Table 1). Furthermore, molecular weights and nitrogen to carbon ratios (N/C) for unique compounds were calculated, and both were overlaid on the van Krevelen diagram (Figures 4C,D). Higher abundance of N-containing compounds (~5% −345 unique compounds out of 6,906 peaks in total) with molecular weight >400 m/z (>400 Da) was observed in bog compared to fen samples. Conversely, fen SOM appeared to have a lower abundance of N-containing compounds (1% −70 unique compounds out of 5,973 peaks in total) and the N-containing compounds in the fen tended to be of lower mass (Supplementary Figure 2).

DBE-O was calculated for all unique fen and bog S- and N-containing compounds to better characterize the OM in the bog and fen samples (Supplementary Figure 3). Interestingly, for S-containing compounds, average DBE-O value for bogs and fens were significantly different (p < 0.01). Average bog DBE-O was 1.36 (±7.03) compared to the fen DBE-O of 4.18 (±8.05), suggesting that S-containing compounds in bog are more oxygenated with potentially higher number of carboxyl groups and fewer C=C double bonds. Average DBE-O values calculated for N-containing compounds were not significantly different between the bog and fen samples.

We used multiple statistical tests to look at changes in SOM primary (central) metabolite composition between bog and fen samples (Figure 5; Table 4) where we observed statistically significant differences in bog and fen primary metabolite abundance. Principal component analysis (Figure 5A) revealed distinct differences in metabolite abundance (explained by PC1, 47.1%) between bog and fen samples. Effect of depth was also explained by PC2 (19.9%). Correlation coefficients of 25 primary metabolite abundances between bog and fen samples (regardless of depth) using Pearson r distance measure were determined at p < 0.05 (Figure 5B). Sugar acids such as D-galacturonic acid and simple sugars such as D-glucose, melibiose, sucrose, fructose, and D-mannose were significantly upregulated (more abundant) in bog samples (regardless of depth) whereas upregulated metabolites in fen samples were mainly classified as acids such as Erythronic acid, glycolic acid, pyruvic acid and amino acids such as glutamine and isoleucine. Another way to view this data is through the use of heatmaps. A heatmap of the same 25 identified metabolites (t-test significance) showed a similar separation (unpregulated sugar acids and simple sugars in bog and organic acid and amino acids in fen) between metabolites identified in different bog and fen samples (Figure 5D). Figure 5C features the top four upregulated metabolites in bog and fen samples created via fold change (FC) analysis and t-test on absolute value changes between metabolites identified in bog and fen samples.

Consistent with changes in SOM and primary metabolite composition, transformation analysis indicated that the biochemical processes associated with the SOM degradation were significantly different between the bog and fen sites (Figure 6A). We observed higher abundances of sugar-associated transformations for the bog sites compared to the fen sites, which generally correlated with carbon transformations (Figure 6A). These data are consistent with Figures 2, 4B where bog sites had a lower abundance of carbohydrate molecules with mass >400 Da but a higher abundance of mono and disaccharides, reflective of decomposition of complex carbohydrates into mono, and disaccharides. However, based on their higher relative abundance in the porewater (Figure 5), these sugars are likely being slowly utilized by microbial communities and/or undergoing slow decomposition. In contrast, fen samples had a higher abundance of complex carbohydrate molecules, but lower abundance of simple, mono, and disaccharides. Despite the higher abundance of complex carbohydrate compounds in the fen site and the higher abundance of sugar - associated transformations in the bog compared to the fen, principal component analysis on specific sugar -associated metabolic transformations on all elemental formulas detected in the bog and fen samples (Figure 6B) indicated similar sugar transformations, suggesting that microbes at both sites have the capability of degrading simple and complex sugar molecules. These data suggest that simple sugars are being utilized at a higher rate in the fen sites compared to the bog sites. The slow decomposition at the bog sites and/or slow uptake of glucose by microbial communities for fermentation or anaerobic cellular respiration suggest additional controls (e.g., thermodynamic, abiotic processes) on SOM decomposition in the bog sites. The high amino acid-associated transformations in the fen samples (Figure 6A) is consistent with a higher abundance of amino acid compounds at the fen sites, reflective of protein degradation and /or release of amino acids by fen vegetation through roots.

Figure 7 and Supplementary Figure 4 show the correlation between network analysis parameters to pH, DOC (mM), and δ¹³C -CH₄ of bog and fen porewater samples. Significant correlations were observed between clustering coefficient and average number of neighbors for SOM molecules to porewater pH, DOC (mM) and δ¹³C -CH₄ values. SOM clustering coefficient increased with increasing pH (Figure 7A), decreasing DOC (Figure 7B), and enriched δ¹³C -CH₄ (Figure 7C) values as in the case of fen and decreased with decreasing pH (Figure 7A) and depleted δ¹³C -CH₄ values (Figure 7C) as in the case of bog. On the other hand, SOM average number of neighbors had an opposite trend to the clustering coefficient (Supplementary Figure 4). Significant correlations (p < 0.05) were also observed between different biogeochemical properties including concentrations of DOC, TN, DIC (CO₂), CH₄, ammonia, nitrate, sulfate, δ¹³Cpw-DIC, δ¹³Cpw-CH₄ (Supplementary Table 4), revealing clear coupling between different biogeochemical properties and GHG emissions.