Microbial Community Analyses Inform Geochemical Reaction Network Models for Predicting Pathways of Greenhouse Gas Production

The mechanisms, pathways, and rates of CO 2 and CH 4 production are central to understanding carbon cycling and greenhouse gas ﬂux in wetlands. Thawing permafrost the terminal electron acceptors are depleted and the system becomes increasingly methanogenic. This suggests that as permafrost regions thaw and dry palsas transition into wet fens, CH 4 emissions will rise, increasing the warming potential of these systems and accelerating climate warming feedbacks.

The mechanisms, pathways, and rates of CO 2 and CH 4 production are central to understanding carbon cycling and greenhouse gas flux in wetlands. Thawing permafrost regions are of particular interest because they are disproportionally affected by climate warming and store large reservoirs of organic C that may be readily converted to CO 2 and CH 4 upon thaw. This conversion is accomplished by a community of microorganisms interacting in complex ways to transform large organic compounds into fatty acids and ultimately CO 2 and CH 4 . While the central role of microbes in this process is well-known, geochemical rate models rarely integrate microbiological information. Herein, we expanded the geochemical rate model of Neumann et al., (2016, Biogeochemistry 127: 57-87) to incorporate a Bayesian probability analysis and applied the result to quantifying rates of CO 2 , CH 4 , and acetate production in closed-system incubations of peat collected from three habitats along a permafrost thaw gradient. The goals of this analysis were twofold. First, we integrated microbial community analyses with geochemical rate modeling by using microbial data to inform the best model choice among equally mathematically feasible model variants. Second, based on model results, we described changes in organic carbon transformation among habitats to understand the changing pathways of greenhouse gas production along the permafrost thaw gradient. We found that acetoclasty, hydrogenotrophy, CO 2 production, and homoacetogenesis were the important reactions in this system, with little evidence for anaerobic CH 4 oxidation. There was a distinct transition in the reactions across the thaw gradient. The collapsed palsa stage presents an initial disequilibrium where the abrupt (physically and temporally) change in elevation introduces freshly fixed carbon into anoxic conditions then fermentation products build up over time as the system transitions through the acid phase and electron acceptors are depleted. In the bog, fermentation slows, while methanogenesis increases. In the fully thawed fen, most of

INTRODUCTION
More than 1650 Pg of organic carbon (OC) is currently sequestered in boreal regions, where it is frozen in permafrost unavailable for microbial decomposition (Hugelius et al., 2013). Climate change-and the resulting warmer temperatureleads to permafrost thaw which makes previously stable OC available to microbial decomposition (Schädel et al., 2016). The decomposition of OC by microorganisms produces the greenhouse gasses CO 2 and CH 4 , which can then escape to the atmosphere exacerbating warming effects in a positive feedback. However, the magnitude of C released and the relative proportion of each gas (which have very different radiative forcings) remains poorly constrained in global models, contributing to uncertainty in climate trajectories (Frolking et al., 2006). Constraining these rates is critical for improving global climate change models and requires understanding of the geochemical and microbiological controls on CO 2 and CH 4 production rates.
Stordalen Mire, Sweden, lies above the Arctic Circle, in the discontinuous permafrost zone, in a region that is disparately affected by climate warming (Kattsov et al., 2005). Indeed, recent studies have highlighted accelerated permafrost thaw within the mire which, over the recent decades, is visually observed as changes in habitat type, hydrology, and vegetation cover (Johansson et al., 2006;Kokfelt et al., 2009). The portions of Stordalen Mire peat that are underlain by intact permafrost are elevated as palsas which are populated by shrubs and woody vegetation. Upon thawing, the elevated palsas collapse and flood creating small thaw pools (Johansson et al., 2006). If these pools remain disconnected from the surrounding water table, the ingrowth of Sphagnum moss over time can elevate the collapse pool into an ombrotrophic Sphagnum bog. In some cases, the insulating effects of the Sphagnum cover can allow permafrost to reform, but further warming eventually causes further subsidence and flooding as the landscape falls to the level of the surrounding water table. Alternatively, thawing of palsa directly adjacent to a fen, or the rapid thawing of palsas surrounding a collapse pool, can connect the collapse feature to the surrounding water table before ingrowth of Sphagnum can occur. All of these processes ultimately lead to increasing minerotrophy as the site transitions into a poor and then a rich fen, leading to a shift in vegetation to sedges and other aquatic macrophytes, characteristic of fens (Zoltai, 1993;Vitt et al., 1994;Jorgenson et al., 2001;Malmer et al., 2005) (Figure 1).
Carbon (C) sequestration in peatlands occurs when net primary production (NPP), measured as CO 2 uptake, exceeds net ecosystem respiration (NER), measured as CO 2 and CH 4 emission, forming mostly cellulosic peat (Clymo, 1987) that accumulates over millennia. In the palsas, dry frozen conditions result in low NPP and NER, such that C sequestration rates are generally low. Sphagnum is thought to produce lower quality litter than either the woody vegetation that precedes it in the thaw gradient or the graminoid vegetation that follows in the fully inundated sites (Malmer et al., 2005) and since low-quality litter contributes to low NER, areas with extensive Sphagnum coverage have the highest C sequestration rates. Fens, with their highly productive yet easily decomposable vascular plant vegetation, have the highest NPP in Stordalen Mire. Thus, across the thaw gradient there are large variations in NER, NPP, and C sequestration rates. Since the 1970s, the areal extent of Sphagnum carpet (combined moss hummocks, semiwet, and non-graminoid wet sites) in Stordalen Mire has declined (Malmer et al., 2005), coincident with overall wetter conditions and an increase in the areal extent of fen habitats (Kokfelt et al., 2009;Bäckstrand et al., 2010).
In addition to differences in NER and the quality of OC inputs, the gradient in habitat types created by thawing across the mire coincides with very different geochemical conditions under which OC is degraded. Cumulatively, these differences result in increasing total CO 2 and CH 4 production rates and decreasing CO 2 :CH 4 ratios across the thaw gradient Hodgkins et al., 2015). In the relatively dry elevated palsas, O 2 easily penetrates the top layer suppressing methanogenesis and oxidizing what little CH 4 is produced before it can be emitted. The net result is that palsas are very weak net C sinks (Johansson et al., 2006) and may sometimes act as net C sources (Bäckstrand et al., 2010) due to interannual changes in precipitation or hydrology that result in wetter or dryer conditions (e.g., Euskirchen et al., 2012). The wetter conditions of the collapsed palsas and bogs allows some methanogenesis to occur, but rates are generally low compared to more fully thawed and inundated fens such that in the collapsed palsa and bogs CO 2 production exceeds CH 4 production . The fully thawed, inundated conditions of the fen stimulate NPP, but also favor methanogenesis with the net result that these habitats emit nearly 1/3 of their net CO 2 uptake as CH 4 (Johansson et al., 2006).
The transition of an ecosystem between a net C sink and a net C source depends on the balance of NER:NPP, however, because of the higher climate forcing potential of CH 4, a habitat that emits CH 4 can exert net positive radiative forcing while still remaining a net C sink. CH 4 has a shorter half-life than CO 2 in the atmosphere and therefore net radiative forcing is timescale-dependent, but for humans that live on an ∼100year timescale, CH 4 forcing remains a relevant consideration. Even on shorter timescales, net radiative forcing could push the ecosystem past a threshold tipping point (e.g., Eppinga et al., 2009) which could have significant impacts including FIGURE 1 | Schematic of the Stordalen chronosequence from the dry palsas sitting atop intact permafrost through recently thawed, inundated, collapse features, more mature, but drier, bogs and fully thawed inundated fens. This figure is from Hodgkins (2016) (used with permission). rapid and irreversible shifts in vegetation and landscape structure. Thus, the CO 2 :CH 4 production ratio is an important consideration for understanding climate feedbacks with these systems (Frolking et al., 2006;Neubauer and Megonigal, 2015;Wilson et al., 2017a). CO 2 and CH 4 production are mediated by complex microbial transformations that, thus far, are still poorly characterized because of their great diversity. For example, homoacetogenesis (i.e., acetate production via the acetyl-CoA pathway: 2CO 2 + 4H 2 → CH 3 COOH + 2H 2 O) (Diekert and Wohlfarth, 1994a,b;Drake et al., 2006) has typically been considered negligible in peatlands (Megonigal et al., 2004;Drake et al., 2006;Tveit et al., 2015) in part because hydrogenotrophic methanogens are thought to outcompete homoacetogens for H 2 substrate under the anaerobic conditions of the peat catotelm (Lovley and Klug, 1983;Jones and Simon, 1985;Hoehler et al., 1999). Recently, however, significant homoacetogenesis rates have been measured in incubations of bog peat suggesting that this pathway of CO 2 consumption may be more important than previously considered (Ye et al., 2014). This may be especially true in cool boreal systems (<10 • C) where the higher growth rates of homoacetogens compared to methanogens may allow them to successfully compete for substrate (Conrad et al., 1989;Conrad and Wetter, 1990;Kotsyurbenko et al., 1993Kotsyurbenko et al., , 1996Kotsyurbenko et al., , 2001Nozhevnikova et al., 1994;Liu and Conrad, 2011). Although, Ye et al. (2014) measured appreciable homoacetogenesis rates via isotopic methods, an alternative explanation of their results is that isotopic exchange between 14 CO 2 and acetate (Eikmanns and Thauer, 1984) may have made it appear as if the microorganisms were incorporating CO 2 into acetate at a greater rate than was actually occurring. In any case, homoacetogenesis only affects the ultimate CO 2 :CH 4 production of the system if the produced acetate is not itself disproportionated into CO 2 and CH 4 , but instead accumulates-which does occur in some habitats . If, on the other hand, the resulting acetate is disproportionated (equation 3), then the net overall reaction is equimolar production of CO 2 and CH 4 : With hydrogenotrophic methanogenesis as the sole H 2 consumption mechanism (i.e., no homoacetogenesis) the set of reactions are: While it is generally thought that hydrogenotrophic methanogenesis will outcompete homoacetogenesis (except possibly at cooler temperatures, Ye et al., 2012), in cases where the processes do coexist, the equations may be written:. In each of these three cases, regardless of the reactions involved, the overall net result is the same: equimolar CO 2 and CH 4 production. However, because these processes fractionate the 13 C in CO 2 and CH 4 to different degrees, the relative importance of these pathways to overall CO 2 and CH 4 production can influence interpretation of isotope-based mass balance models (e.g., Shoemaker and Schrag, 2010;Corbett et al., 2013Corbett et al., , 2015Neumann et al., 2016). Anaerobic decomposition begins with fermentation (which produces CO 2 ) and proceeds through syntrophic reactions and then methanogenesis (which consumes CO 2 while producing CH 4 ). While fermentation does not isotopically fractionate the CO 2 relative to the starting organic material, the net CO 2 produced by the sum of fermentation followed by methanogenesis is fractionated. This fractionation occurs because both methanogenic CO 2 reduction and acetate disproportionation to CO 2 and CH 4 experience a kinetic isotope effect (Whiticar, 1999). Production of 13 C depleted CH 4 leaves the remaining CO 2 pool 13 C enriched relative to the starting organic material (Corbett et al., 2012). Because in a purely methanogenic system, CO 2 and CH 4 are produced in equal proportions (Conrad, 1999), excess non-fractionated (relative to the starting organic material) CO 2 production signifies the presence of alternative electron acceptors and that the system is not at methanogenic steadystate (e.g., continual addition of new bioavailable substrate, biomass production, or a time lag in response to product buildup). This stage is referred to as the acid phase in the waste decomposition literature (O'leary and Tchobanoglous, 2002) due to the high production of volatile fatty acids (VFAs) at this stage of decomposition. As the concentration of VFAs increase in the system, however, thermodynamic equilibrium shifts, increasingly favoring methanogenesis which draws down VFA concentrations until the system reaches a methanogenic steady-state balancing the production of VFAs by the production of CH 4 to maintain favorable H 2 (aq). Differences in CH 4 production pathways among the different thaw habitats in the mire may reflect different thermodynamic controls resulting from the changing geochemical conditions across the sites or perhaps changes in the microbial communities mediating these processes. Understanding the source of variation in reaction rates among the different habitats will be critical for modeling CO 2 and CH 4 production rates from thawing permafrost. To better understand these controls and predict their potential response to further warming, we integrated microbial community analyses with geochemical models to identify the important underlying mechanisms giving rise to CO 2 and CH 4 production dynamics and to calculate the rates and pathways of CO 2 and CH 4 production with transition across the thawing permafrost landscape. To accomplish this, we first characterized the rates, ratios and stable isotopic composition of CO 2 and CH 4 production in incubations of peat collected from habitats undergoing sequential stages of permafrost thaw. In parallel with these gas measurements, we also conducted microbial community profiling of the incubated peat. We then applied a reaction-network model to the CO 2 and CH 4 measurements in these incubations and used changes in the microbial community composition over time to inform our choice of the best model to describe changes in OC decomposition with increasing thaw maturity.

Peat Collection and Incubation
Peat was collected from three sites spanning the thaw gradient at Stordalen Mire, Sweden (68.35 • N, 19.05 • E): collapsed palsa, bog, and fen habitat types. These sites correspond with the sites named "PHB" (collapsed palsa); "S, " "Sphagnum, " or "Bog" (bog); and "E, " "Eriophorum, " or "Fen" (fen) in previous studies at this site (Hodgkins et al., , 2015McCalley et al., 2014;Mondav et al., 2014;Emerson et al., 2018;Singleton et al., 2018;Woodcroft et al., 2018). Peat (320 g) was collected from each of 8-18 cm and 28-38 cm in the collapsed palsa and bog sites, and from 8 to 18 cm in the fen site. It was not possible to recover peat from deeper in the fen site at the time of sampling because of limitations in the coring device. Peat was stored at 2 • C in plastic resealable bags for approximately 4 months before initiation of the incubation experiments. Previous work in our lab demonstrates that incubations of peat treated in this way largely reflect the geochemical and microbiological conditions of the field at the scale of interest in this study and are, therefore, an excellent proxy for studying the controls on CO 2 and CH 4 dynamics under controlled conditions (Hodgkins et al., 2015;Holmes et al., 2015).
Incubations were conducted at 18 • C following previously established protocols Hodgkins et al., 2014Hodgkins et al., , 2015Tfaily et al., 2014;Holmes et al., 2015;Wilson et al., 2017b). Briefly, peat from each site/depth combination was gently homogenized by hand and then peat was divided into eight separate serum vials (40 g peat in each of 125 mL vials) to make replicates for each incubation experiment. 40 mL of degassed, distilled water was added to each incubation vial and then the vials were sealed with a gas-tight pierceable septum. Each vial was flushed with pure N 2 until the headspace measured <0.5% CO 2 via GC-FID. Vials were then pre-incubated in the dark for 25 days to allow reduction of any oxidized species created during the sampling and handling of the peat. After 25 days the vials were again flushed with pure N 2 until the headspace measured <0.5% CO 2 . This procedure has been previously shown to be effective at reducing any oxidized species created during the sampling and handling of the peat that could potentially act as terminal electron acceptors, thus ensuring anaerobic methanogenic conditions (Ye et al., 2012). After the second N 2 flush, the incubation experiments were considered to have started and this time point constituted day 0 of the incubation experiments described in this study. CO 2 and CH 4 concentrations were measured daily by GC-FID coupled to a methanizer for the first 5 days until CH 4 concentrations were high enough to be measured simultaneously with isotopic composition via a Delta V Finnegan isotope ratio mass spectrometer. Then isotopes and concentrations for both CO 2 and CH 4 were collected approximately every 3-5 days for the 50-day duration of the experiment. Since sample sizes for direct injections were 50-250 µL out of the total 45 mL headspace (i.e., <1%) we did not do any gas replacement. CO 2 and CH 4 concentrations were calculated using the respective temperature-corrected Henry's Law constant to account for dissolved as well as headspace components. Replicates were sacrificed on days 0, 25, and 50 and peat from those replicates was stored in LifeGuard (MoBio) at −80 • C for subsequent microbial community analyses.

Isotope Modeling
For the modeling, we built upon the conceptual reaction network rate model of Neumann et al. (2016) to calculate reaction rates in our incubations. The Neumann et al. (2016) model quantifies reaction rates from changes in the concentrations and stable carbon isotopes of CO 2 and CH 4 resulting from the microbial transformations taking place (Shoemaker and Schrag, 2010;Corbett et al., 2012Corbett et al., , 2015Shoemaker et al., 2012;Holmes et al., 2015;Neumann et al., 2016). During C transformation, microorganisms preferentially use the lighter isotope form leaving behind substrate that is increasingly enriched in the heavy isotope. The extent of the isotopic fractionation between substrates and products varies among the different biochemical pathways. There is little to no fractionation between organic C substrates and anaerobic fermentation or respiration products (Lapham et al., 1999), but the fractionation between CO 2 and CH 4 can be as much as 80 depending on the reactions mediating the transformation (Whiticar, 1999). To estimate microbial rates using temporal changes in the stable C isotope ratio (δ 13 C) of CH 4 and CO 2 , we considered five possible anaerobic reactions for OC transformation ( Table 1). Table 1 provides a list of the modeled reactions as well as the fractionation factors associated with the reactions that we used in the modeling. We modeled the rates of acetoclastic methanogenesis, hydrogenotrophic methanogenesis, non-fractionated CO 2 production (i.e., CO 2 resulting from respiration, biomass increase, fermentation, or other non-fractioning mechanism), methane oxidation, and homoacetogenesis in each of the incubations over time and leaving out fugitive CH 4 since the closed incubations did not allow gasses to escape.
Expanding on Neumann et al. (2016), we took a Bayesian approach to specifying parameter values and estimating microbial rates from temporal changes in CO 2 , CH 4 , δ 13 C-CO 2 , Methane oxidation CH 4 + TEA i → CO 2 + reduced TEA i 5-31 e CH 4 consumption * Conceivably the term 'hydrogenotrophy' could also be applied to homoacetogenesis, because that reaction also consumes H 2 , but for the purposes of this paper we use 'hydrogenotrophy' to refer exclusively to the CO 2 reducing pathway of methanogenesis. a Fractionation factor ranges are from Neumann et al. (2016) based on references indicated given in the following notation: ε c = (δ 13 C substrate − δ 13 C product ). b (Gelwicks et al., 1994;Whiticar, 1999;Valentine et al., 2004;Londry et al., 2008). c (Games et al., 1978;Botz et al., 1996;Whiticar, 1999;Valentine et al., 2004;Londry et al., 2008). d (Blaser et al., 2013). e (Whiticar, 1999;Feisthauer et al., 2011). and δ 13 C-CH 4 data from the incubations. More than one million Monte Carlo simulations were performed with varying reaction rates (0-5,000 µmol L −1 day −1 ) and C fractionation factors (ranges given in Table 1) for each time point during each incubation. Initially, the entire listed ranges of reaction rates and fractionation factors were randomly sampled using a uniform distribution across the range. The likelihood of each parameter combination was then calculated by comparing simulated concentrations and isotopic composition to those measured, assuming measurement errors were normally distributed. Likelihood values were used to calculate the joint probability of the simulation describing the data. The probability of each simulation was normalized by the sum of probabilities from all one million simulations, giving each simulation a normalized value between zero and one. Initially, when the entire range of parameter values (rates and fractionation values) were sampled, only a few simulations accounted for over 99% of the normalized probability, indicating that the parameter ranges were too wide to fully probe with one million simulations. In this instance, ranges for rates and fractionation factors were reduced to the minimum and maximum value associated with parameter values in the 1,000 simulations that had the highest probability, regardless of how low the normalized probability value actually was for those simulations. One million additional Monte Carlo simulations were then performed sampling from the more constrained rate and fractionation factor ranges. This process was continued until there were at least 100 simulations that made up 99% of the normalized probability. Results from these ≥100 simulations were used generate probability distributions for the microbial rate values.
We modeled the incubation data using four variants of the rate model (Figure 2), leaving out various combinations of reactions to evaluate the effect these changes to the reaction network would have on (1) the ability of the model to converge (i.e., find a solution) and (2) the final rate estimations. In the first variant, we included all five reactions (acetoclasty, hydrogenotrophy, non-fractionated CO 2 production, CH 4 oxidation, and homoacetogenesis) representing the most complex plausible reaction network. In the second modeling variant, we constrained the CH 4 oxidation rate to be zero and only allowed non-fractionated CO 2 production, hydrogenotrophy, acetoclasty, and homoacetogenesis. In the third modeling variant we constrained the homoacetogenesis rate to zero and only allowed non-fractionated CO 2 production, hydrogenotrophy, acetoclasty, and CH 4 oxidation. In the final modeling variant, we constrained both the CH 4 oxidation and homoacetogenesis rates to zero and only allowed nonfractionated CO 2 production, hydrogenotrophy, and acetoclasty in the model. This final model represents the simplest plausible reaction network (e.g., Corbett et al., 2013).

Microbial Community Analyses
Using 16S rRNA gene amplicon sequencing of the Lifeguard preserved peat samples from each sampling time point, we investigated the relative abundances of microbial community members from each site that are taxonomically affiliated with organisms that are known to perform the modeled microbial reactions. Although 16S alone does not inform microbial function, conservation of functional features among relatives is observed, including those related to these specific modeled functions (Martiny et al., 2015). In particular, methanogenesis as a trait is typically conserved at the order or family level (Martiny et al., 2015), although acetogenesis is not (Drake et al., 2013).

Nucleic Acid Extraction
DNA extractions were performed for replicate peat samples from the field, incubation setup, and days 0, 25, and 50 of the incubations, as described by Mondav et al. (2014) and Woodcroft et al. (2018). Briefly, total nucleic acids were extracted from ∼2 g peat sample using the PowerMax Total Nucleic Acid extraction kit (MoBio), retaining the LifeGuard preservation solution during lysis. DNA was purified by RNaseA digestion, phenolchloroform-isoamyl alcohol purified and ethanol precipitated. Approximately 15 ng DNA from each sample was used as a template in PCR reactions. Libraries for low concentration DNA samples were created using 1 ng of DNA with the Nextera XT DNA Sample Preparation Kit.

Quantitative PCR
To measure microbial cell abundances, based on 16S rRNA gene copies, quantitative polymerase chain reaction (qPCR) was conducted for each community. Each reaction contained 5 µl of 2X SYBR Green PCR Master Mix (Applied Biosystems, Carlsbad, CA, United States), 4 µl of template DNA, and 1 µl of primer mix. Bacterial and archaeal 16S rRNA genes were amplified using the primer set 1406F/1525R (0.4 µM, F -GYACWCACCGCCCGT and R -AAGGAGGTGWTCCARCC). For inhibition control samples, the rpsL primer pair (0.2 µM, F -GTAAAGTATGCCGTGTTCGT and R -AGCCTGCTTACGGTCTTTA) was used to amplify Escherichia coli DH10B only. Each sample and standard were run in triplicate, with three dilutions (1/100, 1/500, and 1/1000) and an inhibition control (1/100 dilution of E. coli DH10B genomic DNA spiked into a 1/100 dilution of the sample). E. coli DH10B genomic DNA dilutions of 10 −2 , 10 −3 , 10 −4 , 10 −5 , and 10 −6 of the 20 ng/µl stock solution were used for the standards. The ViiA7 Real-Time PCR System (Applied Biosystems, Carlsbad, CA, United States) was used for the qPCRs, with the following cycling conditions: 10 min at 95 • C, 40 cycles of [15 s at 95 • C, then 20 s at 55 • C, then 30 s at 72 • C]. To produce a melt curve, a cycle of 2 min at 95 • C was run, and then a final cycle of 15 s at 60 • C. The cycle threshold (Ct) values were recorded and analyzed using ViiA7 v1.2 software, and 16S rRNA gene copy numbers were calculated for each sample. These calculations accounted for the genome size (4,686,137 bp) and 16S rRNA gene copy number (7) of the standard.
Preparation of the 16S library was performed as described, using the workflow outlined by Illumina (#15044223 Rev.B). In the 1st stage, PCR products of ∼466 bp were amplified according to the specified workflow with an alteration in polymerase used to substitute Q5 Hot Start High-Fidelity 2X Master Mix (New England Biolabs) in standard PCR conditions. Resulting PCR amplicons were purified using Agencourt AMPure XP beads (Beckman Coulter). Purified DNA was indexed with unique 8 bp barcodes using the Illumina Nextera XT 384 sample Index Kit A-D (Illumina FC-131-1002) in standard PCR conditions with Q5 Hot Start High-Fidelity 2X Master Mix. Indexed amplicons were pooled together in equimolar concentrations and sequenced on MiSeq Sequencing System (Illumina) using paired end sequencing with V3 300 bp chemistry in the Australian Centre for Ecogenomics according to manufacturer's protocol.
iTag Processing 16S amplicons were processed using the Australian Center for Ecogenomics (ACE) mitag pipeline (AMP). The first 20 bases of all fastq files are then trimmed to remove primer sequence, and quality trimmed to remove poor quality sequence using a sliding window of 4 bases with an average base quality above 15 using the software Trimmomatic. All reads are then hard trimmed to 250 bases, and any with less than 250 bases excluded. Fastq files are then converted to fasta files. Fasta files are processed using QIIME's pick_open_reference_otus.py workflow with default parameters (97% similarity) and taxonomy assignment and alignment features suppressed. The resulting OTU table is filtered to remove any OTU with an abundance of less than 0.005%. Representative OTU sequences are then BLASTed against the reference database (Greengenes version 2013/05 for 16S).
In order to correct the amplicon-based relative abundances for variance in genome copy number of the 16S gene between different lineages, both the OTU data and qPCR data were also processed by CopyRighter (Angly et al., 2014).

Identification and Analysis of Microbial Functional Groups
From the 16S rRNA community profiles, lineages known to perform hydrogenotrophic methanogenesis, acetoclastic methanogenesis, and anaerobic CH 4 oxidation were manually identified through literature review and included those identified at this site by previous molecular characterization of Stordalen Mire microbial communities (Mondav et al., 2017;Woodcroft et al., 2018). The copy-corrected relative abundances of these lineages were calculated within each microbial community. Aerobic CH 4 oxidizing lineages (families Methylocystaceae, Methylococcaceae, Beijerinckiaceae, and Hyphomicrobiaceae) detected at Stordalen Mire by Singleton et al. (2018) were also searched for in the incubation communities. Since both methanogenesis and CH 4 oxidation are fairly taxonomically conserved, inferring function based on 16S rRNA-inferred taxonomy is relatively robust. Acetogens, on the other hand, are phylogenetically diverse, spanning at least 20 genera (Drake et al., 2008). Since homoacetogenesis is not a monophyletic trait, a curated list was developed of demonstrated homoacetogenic lineages that have been isolated, cultured, and experimentally observed to form acetate from CO 2 and H 2 (Supplementary Table S1). Putative homoacetogens in this work are lineages affiliated at either the genus or family level with experimentally verified homoacetogenic lineages. One caveat with this analysis is that some sister taxa diverge in their laboratory ability to perform homoacetogenesis, so inferring function based solely on 16S data is less robust than inferences about CH 4 production and oxidation.

CO 2 and CH 4 Production
CO 2 and CH 4 concentrations increased throughout all incubations indicating the production of these gases over time (Figure 3). CO 2 and CH 4 in both depths from the bog and in the fen were linear over time indicating constant production rates. Concentrations in the collapsed palsa increased exponentially over time indicating that production rates increased during the incubation. Isotopic ratios of the CO 2 and CH 4 varied over time in all incubations, though the variation was not consistent among treatments. In some cases, the stable isotope ratio became more enriched over time and in some cases the isotope ratio became more depleted over time. The shifts in the isotope ratio indicate temporal changes in the reaction rates and the complex nature of the changes make it difficult to estimate which reaction rates are changing without using a reaction network model. For example, if this isotopic shift was purely the result of substrate fractionation due to the closed-system incubations, one could expect that the direction of the isotopic fractionation would be the same for each of the habitats (becoming more enriched with 13 C over time), but since the direction of the fractionation differs (e.g., CH 4 becomes lighter while CO 2 becomes heavier in the fen, but both CO 2 and CH 4 become heavier over time in the collapsed palsa), substrate alteration cannot fully explain this result. Implementation of the reaction network model allows us to calculate the relative rates and understand how the C decomposition dynamics are changing over time and across sites.

Modeling Results
Because incubation data were variable both among replicates and with time, we fit the data with linear or polynomial regressions (as indicated on Figure 3) and used the results of the regressions as input data for the model. Regression types were chosen such that linear regressions were preferred, but when linear regressions explained less than 50% of the variation in the data, we then chose exponential regressions, and if those, too explained less than 50% of the data we then chose polynomial regressions. In two cases, (δ 13 CO 2 in collapsed palsa 28-38 cm and δ 13 CH 4 in bog 28-38cm), not even polynomial regressions explained a satisfactory amount of the variation in the data and in those cases, we used the simple mean of all time points.
All four model variants (Figure 2) converged for all sites/depths. All model variants produced similar rates for all included reactions (Figure 4), which contributed to the difficulty in discerning among them. However, some slight differences in rate results among variants were noted. Inclusion of homoacetogenesis and CH 4 oxidation altered the ratio of acetoclasty to hydrogenotrophic methanogenesis in each site ( Figure 5). In particular, excluding the homoacetogenesis term increased the relative proportion of hydrogenotrophy in most incubations.

Microbial Results
We examined microbial functional guilds to identify the most likely reaction networks giving rise to the observed CO 2 and CH 4 concentrations and isotopes, and to inform model choice. The 16S rRNA copy-corrected relative abundances of microbial community members taxonomically affiliated with hydrogenotrophic methanogens, acetoclastic methanogens, methanotrophs, and homoacetogens were summed in the incubations and compared over time within incubations and across habitats (Figure 6 and Supplementary Table S2). qPCR confirmed that 16S rRNA gene copy numbers -a proxy for cell abundances -did not vary appreciably among habitats, and remained within the same order of magnitude over time in each (data not shown). These observations of fairly equal and stable absolute abundances of microbes in these samples strengthen comparisons of their relative abundances.
In almost all samples, the hydrogenotrophic methanogens were the most abundant functional group. No known lineages of anaerobic methane oxidizers were identified in any of the samples, consistent with their absence in profiles of the microbiota in the field (Mondav et al., 2017;Woodcroft et al., FIGURE 3 | Measured incubation data represented by points and regression equations (dotted lines) used as model input parameters. For the two cases where the isotopic signature of the emitted gas (δ 13 CH 4 in bog deep and δ 13 CO 2 collapse palsa deep) did not change regularly with time [i.e., the regression equations explained less than 50% of the variability in the data (R 2 < 0.5)], we used the average value of the data as the model input parameters. 2018). While OTUs belonging to aerobic methane oxidizers from Methylocystaceae and Methylococcaceae were present in the incubations, at approximately constant or slightly declining relative abundances over time, the highly anoxic conditions of the incubations almost certainly precluded any aerobic methane oxidation activity.
OTUs affiliated with known acetogenic lineages at the genus or family level were present in the incubations in varying abundances (Figure 6). Here, "homoacetogen" refers to an organism whose growth is known to be supported solely by CO 2 + H 2 , creating acetate using the acetyl-CoA pathway (Diekert and Wohlfarth, 1994b). Though homoacetogenesis as a trait is not strongly phylogenetically conserved (Drake et al., 2008), OTUs that were taxonomically affiliated with known homoacetogenic lineages (i.e., isolated and experimentally verified to grow on CO2 and H2; Supplementary Table S2) at the genus or family level were identified as putative homoacetogens, only ∼30% of the OTUs in the dataset could be assigned genustaxonomy, compared to ∼70% at the family level. Putative homoacetogens in the incubations were affiliated with each of the following genera: Clostridium, Moorella, or Treponema, within the families Thermoanaerobaceae or Spirochaetaceae. Overall FIGURE 4 | Rate results from the four variants of the reaction network model. All rates are in µmol month -1 . Stars are medians and boxes vertically outline interquartile ranges of the posterior likelihoods. Box widths are staggered so that overlapping ranges can be visually distinguished. Blue circles with dashed lines indicate abundances of microorganisms thought to perform the relevant reaction (i.e., acetoclasts, hydrogenotrophs, or genus-level homoacetogens) plotted against the secondary right-hand axis. Aerobic methane oxidizers are also included on the CH 4 Ox. graph for comparison. Note that the scale on the fen acetoclasty, hydrogenotrophy, and homoacetogenesis results differ from other habitats due to the higher overall rates in this habitat.
in the incubations, lineages that are taxonomically affiliated with these experimentally verified homoacetogens are present in every sample. At the genus level, these putative homoacetogens increase in relative abundance over time in the collapsed palsa and bog samples (Figure 4). The family level putative homoacetogens are most abundant in the 8-18 cm collapsed palsa samples, where they increase notably over time (Figure 6).

Model Variant 2 and Microbial Differences Across Thaw Gradient
No markers for anaerobic CH 4 oxidizers were identified in any of the samples consistent with other observations at Stordalen Mire . Lineages that matched to aerobic CH 4 oxidizers from Methylocystaceae and Methylococcaceae were present in the incubations, but the 25-day pre-incubation period should have ensured anaerobic conditions (Ye et al., 2012), therefore, it's unlikely that the aerobic methanogens were active. Based on these considerations, we rejected the modeling variants that included CH 4 oxidation (variants 1 and 3). To decide between the two remaining model variants (2 and 4) we conducted a regression analysis of the median rates of acetoclasty, hydrogenotrophy, and homoacetogenesis versus the relevant microbial abundances at the corresponding timepoint over all incubations (Figure 7). Regressions of acetoclasty and hydrogenotrophy vs. the relevant microbial abundances were significant (p < 0.05) for both model variants, however, the standard error of the regressions (S.E.R) were smaller in both cases for variant 2 which indicates that variant 2 provided a more precise fit of the observed data (Figure 7). The regression of homoacetogenesis vs. homoacetogen (genus-level) abundance was not significant (p = 0.3) and since homoacetogenesis is FIGURE 5 | Ratio of modeled rates of acetoclasty:hydrogenotrophic methanogenesis for each of the four model variants. Stars indicate medians and box heights outline interquartile ranges of posterior likelihoods. Box widths are staggered so that overlapping ranges can be visually distinguished. Blue symbols with dashed lines indicate the ratio of acetoclast to hydrogenotroph abundances in the incubations. In the collapsed palsa (28-38 cm) the hydrogenotrophy rate at days 25 and 30 was 0 µmol month -1 (i.e., acetoclasty:hydrogenotrophy was undefined) which prevented us from calculating ratios after day 20 for that incubation.  not modeled by variant 4, we could not compare between the two variants.
Given the better regression results (as indicated by smaller standard error) for acetoclasty and hydrogenotrophy provided by variant 2 (Figure 7) and since microbial analysis revealed evidence of putative homoacetogens in the incubations we rejected variant 4 (which did not include homoacetogenesis). Thus, the microbial data allowed us to narrow our reaction network choices to variant 2 for exploring differences in microbial rates among the various sites and across the thaw gradient.

Collapsed Palsa Results
The variant 2 reaction network estimated that methanogenesis in the collapsed palsa was dominated by hydrogenotrophy (i.e., acetoclasty:hydrogenotrophy < 1 in Figure 5) until the last time point in the 8-18 cm collapsed palsa incubations (Figure 5). In agreement, the microbial methanogen community at this site was dominated by hydrogenotrophs in both shallow and deep samples, to the near exclusion of acetoclasts; this was observed at all time points examined (days 0, 25, and 50) (Figure 4).
Of the three studied habitats, the collapsed palsa had the greatest abundance of putative homoacetogens (Figure 6). Estimated rates of homoacetogenesis were higher in the collapsed palsa than in the other sites over the first 25 days of the incubations; rates dropped sharply after day 25 (Figure 4) as signaled by the inflection point in the measured rates and isotope ratios (Figure 3). At the same time, the estimated ratio of hydrogenotrophy:homoacetogenesis in the collapsed palsa increased over time approaching 1:1 (Figure 5) by day 50 in the 8-18 cm depth. On day 25, the rates of acetoclasty, hydrogenotrophy, and homoacetogenesis all decreased in this incubation (Figure 4), apparently signifying that as less acetate was produced from homoacetogenesis, less acetate was available to be consumed by acetoclasty. The microbial data indicated that this shift was accompanied by an increase in hydrogenotrophic methanogens (Figure 4).
Over the first 20 days of the incubations, the model estimated that non-fractioned CO 2 production was highest in the collapsed palsa incubation compared to incubations of material from the other sites, but this rate declined over time (Figure 4). In the beginning of the collapsed palsa incubations, non-fractionated CO 2 dominated total CO 2 production, but over time, the relative proportion of non-fractionated CO 2 declined leading to an increase in the proportion of methanogenic CO 2 for both depths (Figure 4). Concomitantly, the observed and modeled CO 2 :CH 4 ratio declined over time consistent with the system becoming more methanogenic (Figure 8). The change in modeled CO 2 :CH 4 ratios over time was greatest in the collapsed palsa incubations (>45%) compared to other sites (Figure 8).

Sphagnum-Dominated Bog Results
In the bog, estimated rates of acetoclasty decreased over time (until the last time point in the deep bog), while estimated rates of hydrogenotrophic methanogenesis increased over time (Figure 4). The relative abundance of acetoclasts in the bog did not reflect the trend in acetoclasty rates (Figure 4). A low level of acetoclasts were present in the bog incubations, ranging from <5% of methanogens in the deeper section to 8-21% in shallower bog samples (Figure 6). Thus, in this incubation, acetoclast abundance did not reflect the modeled rates of acetoclasty. For hydrogenotrophy, however, which increased over the course of the incubations there was a corresponding increase in the abundance of hydrogenotrophs (Figure 4). Indicating a good agreement between hydrogenotroph abundance and rates of hydrogenotrophic methanogenesis. The ratio of acetoclasty:hydrogenotrophy in the bog declined over time and approached 1 near the end of the incubations. At the same time, the overall the ratio of acetoclast to hydrogenotroph abundances declined from days 0 to 50 (Figure 5). The deep bog samples had lower proportions of acetoclasts:hydrogenotrophs than the shallow bog samples, though it was not clear that the deep bog had correspondingly lower acetoclasty:hydrogenotrophy rates (Figure 5).
Homoacetogenesis rates were lower in the shallow bog on day 0 than in any of the other sites, but the estimated rate increased over time in both the shallow and deep bog incubations achieving the highest homoacetogenesis rate of any of the sites FIGURE 8 | Ratio of modeled rates of CO 2 and CH 4 production for each of the incubations and model variants. Box widths are staggered as in Figures 5, 6, but interquartile ranges of the CO 2 :CH 4 are tightly constrained and often smaller than the size of the median symbols at this resolution. by day 50 (Figure 4). Consistent with this finding, the relative abundance of genus-level homoacetogens in the bog incubations also increased over time (Figure 4). Paucity of data precluded a statistical analysis of this correlation but there was a visual correlation between homoacetogen abundances and modeled rate distributions (Figure 4).
Non-fractionating CO 2 production was low in the bog incubations compared to the collapsed-palsa incubations and was near zero throughout the deeper bog incubations (Figure 4). But, over the course of the shallow bog incubations, rates of non-fractionating CO 2 production increased and acetoclastic methanogenesis rates decreased such that the proportion of CO 2 produced by non-fractionating mechanisms increased over time (Figure 4). This is the only incubation where there was an increase in non-fractionating CO 2 production. In all other incubations the proportion of CO 2 produced by nonfractionating mechanisms was relatively stable or declined over time (Figure 4). Despite this increase in non-fractionated CO 2 production, the ratio of CO 2 :CH 4 production in incubations with material from both bog depths was relatively stable over time, although it was higher in the shallow bog incubation compared to the deep bog incubation (Figure 8) indicating that the increase in non-fractionating CO 2 production was offset by a decrease in methanogenesis rates.
In the shallow bog archaeal community, the hydrogenotrophic methanogens increased in abundance over time and the acetoclastic methanogens decrease (Figure 4) consistent with the overall modeled decline in acetoclasty:hydrogentrophy over time ( Figure 5). The increasing rates of homoacetogenesis over time in the bog incubations was in agreement with the increasing relative abundance of homoacetogens observed over time in the bog incubations (Figure 4).

Fen Results
The fen incubation had the highest estimated rate of methanogenesis of any of the sites (Figure 4). In the fen, even more than in the bog, acetoclastic methanogenesis dominated over hydrogenotrophic methanogenesis (Figure 6) although hydrogenotrophy increased over time (Figure 5). Concurrent with the rates of acetoclasty, the fen had the highest proportion of acetoclastic methanogens of any of the sites; yet the proportion of acetoclasts:hydrogenotrophs decreased during incubation from >0.20 on day 0 to just 0.12 on day 50 (Figure 5). This decrease in acetoclasts coincided with an overall decline in the estimated ratio of acetoclasty:hydrogenotrophy rates over time (Figure 5).
Estimated rates of homoacetogenesis in the fen incubation were similar to the other sites and remained fairly constant over time (Figure 4). There was no increase in homoacetogen abundance in the fen incubations over time. While the abundance of homoacetogens was highest on day 25, the relative abundance declined over time such that over the 50 days of incubation, homoacetogen abundances stayed constant or declined only slightly (Figure 4). The ratio of non-fractionating CO 2 in the fen incubation, and consequently the CO 2 :CH 4 ratio, was the lowest of any incubation suggesting the fen was the most methanogenic of the sites (Figures 4, 8). The ratio of CO 2 :CH 4 produced in the fen incubations declined slightly over time from 1.25 to approaching the predicted 1:1 ratio consistent with a methanogenic steady-state (Figure 8).

Microbial Results Inform Model Choice
Quantifying rates of greenhouse gas (such as CO 2 and CH 4 ) production is useful for understanding the C budgets of peatlands and the changes induced by global and environmental change. Additionally, correctly accounting the factors that alter isotopic distributions in CO 2 and CH 4 will better constrain modeling efforts based on isotope mass balance (e.g., Corbett et al., 2012;McCalley et al., 2014;Neumann et al., 2016) as well as improve global inverse models (e.g., Houweling et al., 2017). All variants of the reaction rate model we applied to model CO 2 and CH 4 production pathways and rates account for the isotopic fractionation that occurs as microorganisms process organic material into CO 2 and CH 4 . Since all four model variants converged (i.e., a solution was found for all 4 model variants), all of the model variants were mathematically plausible and all five pathways (acetoclasty, hydrogenotrophy, homoacetogenesis, CH 4 oxidation, and non-fractionating CO 2 production) are plausibly utilized in wetland systems to some degree. Given the similarity in rate results among the four models, and since it was not possible to mathematically choose the best variant to model these data, we used the microbial data to help inform the best model variant for inferring whether homoacetogenesis and CH 4 oxidation were occurring among the different habitats.
Mathematically, under methanogenic conditions, homoacetogenesis should not ultimately affect net radiative forcing because homoacetogenesis increases the production of acetate which should, in turn, stimulate the relative proportion of the acetoclastic pathway relative to the hydrogenotrophic CH 4 production pathway. However, the net ratio of CO 2 :CH 4 production should be the same whether homoacetogenesis occurs or not (see equations 1-4). In our models we observed this most strongly in the collapsed palsa (both depths) and the bog 8-18 cm incubations where model variants that included homoacetogenesis (variants 1 and 2) had higher acetoclasty:hydrogenotrophy rates than model variants that excluded homoacetogenesis (Figure 5). In the incubations of the deeper bog and fen peat, differences among the model variants were much less pronounced (Figure 5) while the CO 2 :CH 4 ratios were similar among all model variants except in the case of the collapsed palsa where the variant that included both homoacetogenesis and CH 4 oxidation had the highest CO 2 :CH 4 ratio after 25 days of incubation (Figure 8). Given the temperature sensitivity of homoacetogenesis (Ye et al., 2012), it is possible that homoacetogenesis may be disproportionally sensitive to temperature potentially leading to mechanistic shifts in CO 2 and CH 4 production with climate warming that may alter the ecosystem net radiative forcing. Thus, understanding the importance of homoacetogenesis as well as the subsequent fermentation of acetate produced in this way is a critical gap in understanding CO 2 and CH 4 production dynamics. This is especially true in peatlands because more fully thawed habitats (e.g., fens) tend to utilize more acetate in the production of CH 4 than the collapsed palsa and bogs. This trend has been observed at this and other peatlands spanning different continents (Corbett et al., 2012;Tfaily et al., 2013;Hodgkins et al., 2014) and is exemplified in our incubations by the lower acetoclasty:hydrogenotrophy ratios in the bogs and collapsed palsa relative to the fen (Figure 5).
Rates of CH 4 production in the incubations increased along the thaw gradient consistent with observations from field CH 4 fluxes (e.g., McCalley et al., 2014). There was also an increase in total methanogen abundance from palsa to bog to fen (Figure 6). As expected, we found microbial evidence of both hydrogenotrophic and acetoclastic methanogens. There are seven known orders of methanogens within the phylum Euryarchaeota, of which, five were identified in our samples: Methanomicrobiales, Methanobacteriales, Methanocellales, Methanosarcinales, and Methanomassiliicoccales (Vanwonterghem et al., 2016) (Figure 3). The first three are only known to produce CH 4 via the hydrogenotrophic pathway. The fourth order, Methanosarcinales, is acetoclastic, however, only the genus Methanosaeta is obligately acetoclastic while Methanosarcina can perform both hydrogenotrophic and acetoclastic methanogenesis under favorable conditions. The fifth order, Methanomassiliicoccales, is a methylotrophic methanogen that has been observed at low abundance in the field at the fen site (Mondav et al., 2017;Woodcroft et al., 2018). The methanogen community was dominated by hydrogenotrophs in all incubations, with the collapsed palsa being almost exclusively hydrogenotrophic and an increasing presence of acetoclastic methanogens across the thaw gradient. The higher abundance of acetoclasts in the fen environment is consistent with observations at other peatlands that fens, with their more microbially available organic material generally have higher rates of acetoclastic methanogenesis than adjacent bogs (e.g., Chanton et al., 2008).
The presence of aerobic CH 4 oxidizers in the community profiles likely represents inactive community members awaiting reinstatement of oxic condition such as might occur in the field through water table fluctuations or intrusion of sedge roots. Although not conclusive, we have also failed to find an appropriate inorganic terminal electron acceptor to support CH 4 oxidation if it is occurring. Organic matter itself serves as an electron acceptor for microbial respiration in a variety of aquatic environments including peatlands, but the exact mechanism or extent to which this occurs is poorly understood (Lovley et al., 1996;Coates et al., 2002;Kappler et al., 2004;Keller et al., 2009;Roden et al., 2010;Smemo and Yavitt, 2011;Keller and Takagi, 2013;Tfaily et al., 2013Tfaily et al., , 2014Wilson et al., 2017b). Porewater sulfate is below detection while nitrate (Hodgkins pers comm) and iron (Chanton unpublished data) concentrations are typically sub-millimolar to below detection limits in the field at these sites. Additionally, the pre-incubation period in our experiment is sufficient to reduce any small amounts of inorganic TEAs that may have been present before the start of the incubations (Ye et al., 2012). This coupled with the dearth of ANME-1,-2,-3 and Methylomirabilis is inconsistent with a significant contribution of anaerobic methane oxidation in our experiments. Based on these considerations, we eliminated model variants that included CH 4 oxidation.
We did find evidence of putative homoacetogens in the microbial analyses at both the genus and family level. While homoacetogenesis is not well-conserved and therefore inferences about homoacetogenic activity based purely on taxonomic similarity are tenuous, we generally observed good agreement between changes in genus-level homoacetogen abundances and variability in rate results in all incubations (Figure 4) suggesting that in this system, taxonomic similarity usefully identified putative homoacetogens.
Based on the microbial data indicating that homoacetogenesis could occur in the samples, but that anaerobic CH 4 oxidation was unlikely, we chose variant 2, which included all reactions except CH 4 oxidation, as the best variant to describe differences in rates across the thaw gradient. Thus, the microbial data allowed us to eliminate potential reactions from the modeling network when geochemical data could not distinguish among the different variants.

Collapsed Palsa
Peat in permafrost settings is elevated into palsas above the surrounding landscape by ice expansion and accretion. When the supporting ice thaws and drains, the palsa collapses and frequently floods creating a heavily inundated thaw feature. This collapse results in an abrupt transition from elevated aerobic conditions in the surface of the palsa, to depressed inundated, anaerobic conditions in the resulting surface peat. Thus, permafrost collapse features may best be understood in the context of other early stage organic-rich organic environments (such as landfills) (Tchobanoglous and Kreith, 2002). In our collapsed palsa site, acetate concentrations were the highest of any of the thaw features ( Table 2). This may reflect an early acid-phase as is seen in other early stage anaerobic environments such as covered landfills or lake sediments (Burdige, 2006). The identification of these processes as "acid phase" is supported by the increasing rates of both acetoclastic and hydrogenotrophic methanogenesis over time in the collapsed palsa incubations after day 20 (Figure 4) and the decline in CO 2 :CH 4 production rates (Figure 8) suggesting that the system was becoming increasingly methanogenic as the incubations progressed, possibly due to the depletion of microbially available organic polysaccharides and the resulting build-up of acetate. Concomitantly, excess non-fractionated CO 2 (e.g., heterotrophic respiration products) decreased over time suggesting declining fermentation rates.
Initially, homoacetogenesis rates in the collapsed palsa were greater than hydrogenotrophic methanogenesis rates, but this ratio declined over time, until the last time point when hydrogenotrophic methanogenesis rates were greater than homoacetogenesis in the shallow collapsed palsa (Figure 4). This finding was similar to results observed by Ye et al. (2014) for other peatlands and was interpreted to indicate direct competition for H 2 between homoacetogens and hydrogenotrophs. Our results differ from Ye et al. (2014) in that, while hydrogenotrophy became increasingly important relative to homoacetogenesis over time, the absolute rates of both processes also increased over time after day 20. The high homoacetogenesis rates in the collapsed palsa incubation may suggest that the partial pressure of H 2 was elevated (Ye et al., 2014) allowing homoacetogens to effectively compete with hydrogenotrophic methanogens for substrate. If this dynamic is reflective of that occurring in the field at the collapsed palsa site, such elevated H 2 could contribute to the high levels of acetate found at this habitat ( Table 2), and also suggests that something may be inhibiting acetoclastic methanogenesis which was very low compared to the other habitat types (Figure 4). Ye et al. (2014) suggested that declining fermentation rates resulted in a decrease in both acetate and H 2 production which limited both homoacetogenesis and hydrogenotrophic methanogenesis. The increase in the rates of both reactions as well as the increase in homoacetogens throughout our incubations, however, cannot be easily reconciled with this interpretation. Instead, we infer, that the declining excess of non-fractionated CO 2 over time in our incubations supported increasing rates of hydrogenotrophic methanogenesis, which also drew down excess H 2 concentrations. At some point (i.e., approximately 35 days) the H 2 concentrations were low enough that hydrogenotrophy outcompeted homoacetogenesis rates and after that shift, the peats became increasingly methanogenic.
Microbial abundances and productions rates do not always correlate strongly (e.g., collapsed palsa 28-38 cm Figure 4). This may reflect that the microbial community is not comprehensively and exhaustively characterized. For example, large portions of the microbial community have not been isolated in pure culture in the laboratory and their function is unknown, or microorganisms which are classified as homoacetogens may also be able to generate methane acetoclastically even if this activity has not been documented under laboratory culture conditions. Additionally, some microorganisms may alter their activity depending on environmental conditions such as homoacetogens which reverse the reaction to produce CO 2 under certain thereby actually decreasing acetogenesis rates (Pierce et al., 2008). Taken on their own, neither the geochemical data nor the microbial data is capable of fully explaining the system, when considered as complementary data, however, a clearer picture of CO 2 and CH 4 production dynamics emerges.

Sphagnum-Dominated Bog
At Stordalen Mire, the intermediate stage of peatland permafrost thaw is characterized by a Sphagnum-dominated bog feature. This site is elevated relative to the collapsed palsa but still below intact palsa such that the surface of the peat is intermittently inundated. The surface peat at these sites is influenced by the oscillating water table that alters O 2 availability over time, and the release of sphagnum acid from the moss that reduces the pH and inhibits microbial activity. In contrast to the collapsed palsa, the estimated rates of acetoclastic methanogenesis in the bog were higher than that for hydrogenotrophic methanogenesis at all time points. This result is in contrast to what was estimated in the field from isotopic mass balance of chamber flux measurements (McCalley Values are averages ± one standard deviation.  Figure 5). However, the ratio of acetoclasty to hydrogenotrophy in the bog incubation was always lower than in the fen incubation, consistent with field findings that hydrogenotrophy was relatively more important in the bog than in the fen ( Figure 5). However, in our study, the estimated acetoclasty to hydrogenotrophy ratio in the bog incubation declined over time as did the ratio of acetoclastic to hydrogenotrophic methanogens suggesting that the bog incubations were becoming more hydrogenotrophic over time. The shift from acetoclastic to more hydrogenotrophic methanogenesis is consistent with thermodynamic considerations that acetoclastic methanogenesis is more energy yielding than hydrogenotrophic methanogenesis (at standard conditions, Conrad, 1999) thus rates of hydrogenotrophy increase after more energetic pathways have been exhausted. Measured acetate levels in the bog are low: 5-9 µM ( Table 2), which is at the lower limit of what is considered the threshold for acetate uptake even for the acetoclastic methanogens in Methanosaeta (Griffin et al., 1998), thus we suggest that as the incubations progressed, the organic matter became more decomposed, acetoclasty became substrate limited and diagenesis proceeded through the next energetically favorable reaction (i.e., CO 2 reduction). In addition to substrate limitation, the declining ratios of acetoclasty:hydrogenotrophy and corresponding declines in acetoclast:hydrogenotroph abundance in the bog (Figures 5, 6) could indicate inhibition of acetoclastic methanogens in the bog. For example, the acetoclastic methanogen genus Methanosaeta (which was least abundant in the bog) is inhibited by ammonia, high concentrations of VFAs (Karakashev et al., 2006), and glycine (Sherwood, 1972), the latter thought to occur by the competitive uptake of Ca 2+ (Sherwood, 1972). Ammonia concentrations are higher in the bog than in the fen, while Ca 2+ concentrations are lower ( Table 2), possibly indicating inhibition of Methanosaeta in the bog relative to the other sites.

Fen
Fens are typically considered the most mature thaw features at Stordalen Mire . They are characterized by high water table and the predominance of sedge vegetation. Sedges are important conduits for conducting gasses between anaerobic peat depths and the atmosphere. They can transport O 2 to otherwise oxygen-depleted depths possibly promoting CH 4 oxidation and allowing CH 4 and CO 2 to transit directly to the atmosphere (e.g., Chanton, 1991;Whiting and Chanton, 1993;Bhullar et al., 2013). The closed incubations that we used in our study, however, prevent CH 4 loss and so the influence of vegetation as a conduit was controlled for in this experiment. However, differences in the microbial availability of organic matter deposited by sedge vegetation versus Sphagnum mosses still remains. Estimated rates of CH 4 production were highest in the fen incubations compared to the other habitat types, consistent with previous incubation studies (Hodgkins et al., , 2015 and field observations of habitat emissions . While the modeled rates of hydrogenotrophic methanogenesis in the fen were higher than the other sites, the fen incubations also had much higher rates of acetoclasty than the other sites (3-4 times higher during some phases of the incubations), such that the acetoclasty:hydrogenotrophy ratio was highest in the fen of all the sites (Figure 5). This is consistent with observations from other peatlands that the more biologically available organic matter in fens stimulates acetoclastic methanogenesis .
Metagenomic analyses of the microbial carbon degradation capacity in the field at this same site, has indicated that not only does the fen have 2-3 times more microbial biomass per g of peat relative to the other habitats, but that the fen also has the most diverse range of microorganisms capable of degrading polysaccharides (the first step in cellulose decomposition) . While there are still more hydrogenotrophs than acetoclasts in the fen, the fen does have the highest abundance of acetoclasts of all the sites (Figure 6). One possible explanation for the apparent discrepancy between the microbial abundances and the modeled rates is that all microbial acetoclasts and hydrogenotrophs may not yet be comprehensively and exhaustively characterized. There may exist microorganisms accomplishing either of these reactions that have not yet been discovered or microorganisms that have been characterized as hydrogenotrophic may also be able to produce methane acetoclastically, but that function has not yet been documented. However, there have been previous extensive metagenomic studies conducted at our field site from which metagenome-assembled genomes for methanogens, including deeply branching lineages were assembled (Supplementary Table S2; Mondav et al., 2014;Woodcroft et al., 2018). Given this extensive characterization, it is unlikely that there are large numbers of previously undescribed acetoclasts at the site. Alternatively, the mismatch between microbial abundances and reaction rates may have a thermodynamic explanation.
Under the non-standard conditions of our experiment, it is possible that high CO 2 , H 2 and or low acetate concentrations (e.g., see bog and fen concentrations Table 2) may create conditions under which hydrogentrophic methanogenesis actually becomes more energetically favorable than acetoclasty. While, at standard conditions, acetoclasty is a higher energy yielding reaction than hydrogenotrophy (on a mol per mol basis), at high CO 2 , H 2 , and/or low acetate concentrations (e.g., see bog and fen Table 2) hydrogenotrophy can become more energy yielding thereby favoring the proliferation of hydrogenotrophic methanogens and forcing acetoclasts to produce CH 4 at higher rates to fulfill the same energy demands. Additionally, methanogens may simply be poorly characterized and some putative hydrogenotrophs may actually be capable of disproportionating acetate as well. For example, Methanosarcina are capable of acetoclastic methanogenesis, but they can also produce CH 4 via hydrogenotrophy. On the other hand, Methanosaeta are considered obligate acetoclasts (Karakashev et al., 2006) and their exclusive presence in the fen is consistent with the high rates of acetoclasty estimated for the fen incubation relative to the others incubations (Figure 4).
Although hydrogenotrophic CH 4 production was once observed in a culture of Methanosaeta via direct interspecies electron transfer from Geobacter metallireducens to Methanosaeta harundinacea, which allowed the reduction of CO 2 (Rotaru et al., 2014), the Rotaru et al. (2014) study is the only observation of such a process, and in our incubations, there was a near-zero abundance of Geobacter spp. Methanosaeta have lower threshold acetate requirements than Methanosarcina (5-20 µM vs. 1 mM) and are only known to outcompete Methanosarcina at low acetate concentrations (Griffin et al., 1998). Thus, the low acetate concentrations in the fen (11.0 ± 2.8 µM) may have allowed Methanosaeta to compete under the conditions of the fen habitat. The estimated ratio of hydrogenotrophy to homoacetogenesis was consistently high in the fen incubation suggesting that H 2 concentrations were low. Excess unfractionated CO 2 concentrations were also low, and CO 2 :CH 4 production ratios approached 1, suggesting that the system was tending toward a methanogenic steady-state. In this way, the fen most closely approximates the theoretical thermodynamic limitation for TEA-depleted anaerobic cellulose degradation. Cumulatively our results suggest that the organic matter in the fen was the most microbially available of the three sites.

Changes Across the Thaw Gradient
There is evidence in the incubations that the collapsed palsa is not at a theoretical methanogenic steady-state (i.e., CO 2 :CH 4 did not = 1). The high CO 2 :CH 4 ratio in this habitat (Figure 8) suggests that the rate of fermentation is in excess of methanogenesis. High rates of fermentation without subsequent methanogenesis would lead to increased H 2(aq) which in turn would favor homoacetogens (Weimer, 1998;Drake et al., 2006;Ye et al., 2014) explaining the high rates of homoacetogenesis in this habitat type. However, the rates shifted over time such that by the end of the incubations, hydrogenotrophic methanogenesis started to dominate over homoacetogenesis ( Figure 5). These results suggest that initially there were organic electron acceptors present in the collapsed palsa incubations that got used up over the course of the incubations. As respiration rates declined, the rates of both methanogenesis pathways gradually increase in the incubations. This trend is consistent with an exhaustion (titration) of organic electron acceptors and a switch to a more methanogenic system. Without a mechanism for re-oxidation, eventually all terminal electron acceptors in a system will become exhausted-except for CO 2 which is produced during anaerobic fermentation. Thus, any reaction pathway in these closed, anaerobic incubations other than CO 2 reduction is ultimately unsustainable and the system will tend toward methanogenesis. Once all of the other terminal electron acceptors have been reduced, methanogenesis rates must balance fermentation, because fermentation is self-limiting, i.e., the build-up of fermentation products thermodynamically inhibits further fermentation from occurring. In this way, the incubations of the collapsed palsa were reflecting the transition from an acid-phase-dominated (production of electron-rich hydrocarbons, e.g., butyrate, propionate, etc.) environment into the early-stage methanogenic conditions that will dominate the subsequent thaw features.
In the shallow bog incubation, there was always more CO 2 produced than could be accounted for by methanogenesis alone (i.e., CO 2 :CH 4 > 1 Figure 8) suggesting the occurrence of alternative electron acceptors or that the bog, also, was not yet at steady-state (i.e., still in the acid phase). Over the course of the incubations, acetoclastic methanogenesis (particularly at 8-18 cm) decreased slightly corresponding to an increase in the non-fractionated CO 2 production over time (Figure 4). Together, these results suggest that fermentation in the shallow bog incubation was slowed relative to the fen which had already progressed through the fermentation-dominated acid phase and was approaching a methanogenic steady-state. The slowing of fermentation that held the bog in the acid phase possibly resulted from the action of Sphagnum-derived acids, which are thought to inhibit the microorganisms responsible for breaking down complex cellulose into the substrates needed for fermentation .
In the fen incubation, the estimated CO 2 :CH 4 ratio declined slightly but it was relatively stable near the predicted 1:1 ratio consistent with the fen being at methanogenic steady-state. The fen feature is similar to the collapsed palsa feature in that it is fully inundated and it has an abundance of sedge vegetation, fens develop as the active layer deepens over time (more thawing) and eventually the feature floods becoming fully inundated and connected with the surrounding water table. Thus, the fen is a more mature thaw feature than the collapsed palsa. Over time the collapsed palsa incubations proceed through the acid phase and toward the end of the incubations, approach methanogenic steady state, whereas the more mature fen material had already achieved this steady state by the start of the incubations. We suggest this progression through the acid-phase as a possible conceptual model for the field thaw transition: the newly thawed immature collapse features have high fermentation and homoacetogenesis rates as fresh organic material is suddenly subjected to anaerobic conditions. The excess substrate is rapidly depleted, however, and eventually, in the absence of other TEAs or further inputs, the system progresses toward a methanogenic steady-state as the thaw site matures. This would suggest that upon thawing initially there may be a rapid release of CO 2 coincident with the high fermentation rates, but over time (decadal timescale) there will be a shift toward increasingly methanogenic conditions. Because CO 2 and CH 4 have very different global warming potentials, the implications of this is that models may need to account for maturity of the thaw feature when estimating climate forcing.

How Does Model Variant Choice Influence the Conclusions?
The differences in the rates among the various modeling variants were subtle and overall trends across the thaw gradient discussed above remain unchanged regardless of the model variant choice. The collapsed palsa had the highest estimated CO 2 :CH 4 ratio and the fen had the lowest estimate CO 2 :CH 4 ratio, regardless of the model variant chosen. The relative importance of acetoclasty increased across the thaw gradient regardless of model variant. Finally, CO 2 :CH 4 decreased across the thaw gradient in all variant conditions. Thus, the overall conclusions regarding the stages of anaerobic decomposition and relative bioavailability of organic matter are ultimately unaffected by the model variant choice. However, including homoacetogenesis decreased the relative importance of hydrogenotrophy compared to acetoclasty and somewhat increased the excess unfractionated CO 2 production. This change has implications for estimating CO 2 and CH 4 production dynamics from isotope mass balance approaches in the field where, unlike in a closed incubation bottles, methane can easily escape from porewater and complicate rate estimations. In Neumann et al. (2016), inclusion of homoacetogenesis shifted the estimated acetoclastic:hydrogenotrophy methanogenesis pathway ratio from ∼1 to ∼3 for a thermokarst bog in Interior Alaska when the models allowed for 'fugitive methane' or the loss of methane from porewater. Excluding homoacetogenesis from field-based microbial reaction network models results in an increase in the estimation of hydrogenotrophy relative to acetoclastic methanogenesis.
Our modeling results suggested that acetoclasty dominated in the early period of the incubations in particular in the modeling variants that included homoacetogenesis. This is in contrast to previous reports that hydrogenotrophic methanogenesis dominates in situ methane production in the bog , however, in situ measurements did not consider the fractionating effect of homoacetogenesis when using carbon stable isotopes to infer pathways. Similar results to McCalley et al. (2014) suggesting the dominance of hydrogenotrophy in bogs abound outside of Stordalen as well (e.g., Chanton et al., 1995;Chasar et al., 2000;Duddleston et al., 2002;Galand et al., 2005;Keller and Bridgham, 2007;Corbett et al., 2012;Zalman et al., 2018). The trade-off between homoacetogenesis and hydrogenotrophy in our model suggests that previous efforts could have overestimated the importance of hydrogenotrophic methanogenesis through an inadequate accounting of the isotopic fractionation introduced by homoacetogenesis. Understanding all of the microbial mechanisms contributing to isotopic fractionation is important for modeling efforts aimed at understanding environmental change, since each mechanism may be accomplished by different microbial consortia that may respond differently to increasing temperature, drought, or changing oxidation conditions.

CONCLUSION
Estimating microbial rates based purely on CO 2 and CH 4 concentrations and isotopes could not alone distinguish between multiple reaction network variants. Neumann et al. (2016) similarly could not distinguish among model variants when calculating in situ CO 2 and CH 4 production rates in Alaska where the fugitive CH 4 term was unknown, but assumed to be large. Using closed incubations, where there is no fugitive CH 4 term, we were hoping to better constrain the reaction network. While there were small differences in the rates of modeled reactions among variants, these differences were not enough to identify a single best variant. The result is that all reaction networks are equally mathematically plausible. However, by coupling the modeling results with an analysis of the microbial community, we were able to narrow the choices. Because the microbial community suggested homoacetogens were abundant, while CH 4 oxidizers were absent, we were able to identify the variant including homoacetogenesis and excluding CH 4 oxidation as the most plausible for modeling our incubations. We then used the results from this variant to describe differences in OC decomposition across the thaw gradient. Our coupled approach revealed a progression across the thaw gradient both in terms of the stage of anaerobic decomposition as well as the microbial availability of the organic matter. The earliest thaw stage, the collapsed palsa, had typical early acid-phase anaerobic decomposition with the production of excess CO 2 and acetate coupled with high rates of homoacetogenesis. Conditions in the bog appeared to be further along the methanogenic diagenetic gradient but were limited by the biological availability of the available organic matter. Thus, there was a shift over time from acetoclastic toward more hydrogenotrophic methanogenesis. Excess unfractionated CO 2 production (in the shallow bog) was higher than expected for a purely methanogenic system consistent with hydrogenation of the aromatic acids released from Sphagnum which would lower H 2(aq) without producing CH 4 (e.g., Wilson et al., 2017b). The fen had the highest rates of acetoclasty and a nearly 1:1 ratio of CO 2 :CH 4 production suggesting that the fen had the most microbially available organic matter and was nearest a methanogenic steady-state.
As permafrost peatlands thaw, changes in both the dominant vegetation and the geochemical conditions result in changing CO 2 and CH 4 production dynamics. Cumulatively, our results suggest that upon permafrost thawing there may be an initial rapid release of CO 2 , but that as the system transitions toward increasingly warm wet habitat, there is a concomitant shift from CO 2 to CH 4 production. Because of the different global warming potentials of these two products (CO 2 and CH 4 ) the implication of this finding is that models may need to account for maturity of the thaw feature when estimating climate forcing.
Homoacetogenesis, which is typically ignored in peatlands, appears from both the microbiology and the rate modeling results to be a critical pathway for CO 2 consumption. Ignoring this reaction, which alters the isotopic signature of the CO 2 (Neumann et al., 2016), leads to overestimation of the importance of hydrogenotrophic methanogenesis by isotopic mass balance models. Because hydrogenotrophy is the least energetic pathway in anaerobic C degradation, it may be the most sensitive to temperature variation, thus overestimating this pathway may lead to erroneous conclusions regarding the sensitivity of the system to climate warming. Alternatively, warmer temperatures may shift the thermodynamics of the system to favor hydrogenotrophy over homoacetogenesis (Ye et al., 2014) exacerbating the warming effects of CH 4 production. Additional work applying this reaction network model in the field and integrating it further with microbial analyses will be key to fully elucidating the critical C cycling pathways and understanding the feedbacks between thawing permafrost landscapes and the atmosphere.

DATA AVAILABILITY
All data is publicly available via the IsoGenie database (https: //isogenie-db.asc.ohio-state.edu/). The raw data and MATLAB code for the model supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
RW and JC conceived the experiments, contributed to model development, analyzed the data, and drafted the manuscript. RN conceived the experiments, developed and interpreted model results, and contributed to manuscript writing. BW, KC, and NR analyzed microbial data and contributed to writing the manuscript. SH conducted geochemical analyses and contributed to writing the manuscript. SS conceived the experiments and contributed to the manuscript. BB and BW analyzed and interpreted microbial data. GT and VR conceived the experiments, analyzed and interpreted microbial data, and contributed to the manuscript.