Opposing Effects of Climate and Permafrost Thaw on CH4 and CO2 Emissions From Northern Lakes

Small, organic‐rich lakes are important sources of methane (CH4) and carbon dioxide (CO2) to the atmosphere, yet the sensitivity of emissions to climate warming is poorly constrained and potentially influenced by permafrost thaw. Here, we monitored emissions from 20 peatland lakes across a 1,600 km permafrost transect in boreal western Canada. Contrary to expectations, we observed a shift from source to sink of CO2 for lakes warmer regions, driven by greater primary productivity associated with greater hydrological connectivity to lakes and nutrient availability in the absence of permafrost. Conversely, an 8‐fold increase in CH4 emissions in warmer regions was associated with water temperature and shifts in microbial communities and dominant anaerobic processes. Our results suggest that the net radiative forcing from altered greenhouse gas emissions of northern peatland lakes this century will be dominated by increasing CH4 emissions and only partially offset by reduced CO2 emissions.

2 of 16 CH 4 emissions (Yvon-Durocher et al., 2011. Temperature sensitivity of net CH 4 emissions is however often assessed within individual ecosystems, for example, seasonally, and it is uncertain whether this sensitivity reflects spatial differences among lakes at broad geographical scales (Morrissey & Livingston, 1992;Rasilo et al., 2015;Sepulveda-Jauregui et al., 2015). Heterotrophic respiration and net emission of CO 2 are also temperature-dependent (Gudasz et al., 2010;Kosten et al., 2010), but relationships between lake net CO 2 emissions and temperatures are often less clear than for CH 4 (Lapierre et al., 2015;Sobek et al., 2005). While most boreal peatland lakes are net sources of CO 2 (Hastie et al., 2018;Lapierre et al., 2015;Weyhenmeyer et al., 2015), CO 2 emissions are influenced not only by respiration but also by CO 2 uptake through primary productivity, which in turn is influenced by nutrient availability and water column light conditions (Lapierre & Del Giorgio, 2012;Pacheco et al., 2014).
Permafrost thaw in peatland regions may have additional direct and indirect effects on lake CO 2 and CH 4 emissions through several processes. Permafrost thaw along lake edges may directly impact GHG emissions through abrupt thaw and the deposition of organic matter into the lakes or through the gentle expansion of saturated lake areas (Turetsky et al., 2020;Walter Anthony et al., 2018). More indirectly, thaw opens both shallow and deep groundwater connections, bringing water to lakes with high concentrations of dissolved organic matter (DOC), nutrients, and weathering products, respectively (Frey & McClelland, 2009;Laudon et al., 2012). Within-lake primary productivity increases with greater nutrient availability, often associated with groundwater inputs but can be inhibited if concentrations of DOM cause substantial light attenuation (Bogard et al., 2020). Greater delivery of CO 2 and CH 4 produced in surrounding organic soils can also cause high lake emissions (Paytan et al., 2015;Weyhenmeyer et al., 2015), while connectivity to deeper groundwater may increase pH and reduce lake CO 2 emissions through pH-driven changes in carbonate equilibrium (Finlay et al., 2015). Higher CH 4 emissions may also be driven by the availability of labile organic matter, supplied from the surrounding landscape, littoral emergent vegetation, or from within-lake primary productivity (Bastviken et al., 2004;DelSontro et al., 2016;Emilson et al., 2018). Conversely, increased delivery of terminal electron acceptors, for example, sulfate and iron, may inhibit methanogenesis (Sivan et al., 2011;Yu et al., 2016). Indirect effects of climate warming through permafrost thaw may thus either enhance or dampen the direct effects of climate warming leading to unclear climate feedbacks from a broader scope than just temperature effects.
We measured CO 2 and CH 4 emissions, alongside water chemistry and microbial community dynamics, during two ice-free seasons from 20 shallow lakes across a 1,600 km permafrost transect within the Interior Plains of western Canada (Figure 1). The Interior Plains of western Canada represents a unique opportunity to use a space-for-time approach to assess climate sensitivity of lake CO 2 and CH 4 emissions, as the region has broadly similar vegetation and a similar glacial history throughout the whole transect. The region is characterized by generally flat topography, sedimentary bedrock overlain by thick, heterogeneous surficial geology (Geological Society of Canada, 2014), and widespread peatlands with abundant lakes (Figure 1 insets). The region has a relatively dry climate throughout, with mean annual precipitation between 300 and 450 mm and mean annual runoff between 50 and 150 mm, but it has a large gradient in mean annual average temperatures (MAAT) between −8°C and +2°C (Table S1 in Supporting Information S1). Peatlands in the Interior Plains started developing shortly after deglaciation ∼10,000 years ago, and today have 2-5 m deep peat deposits (Halsey et al., 1998;Hugelius et al., 2020). Permafrost expanded southward after the Holocene thermal maximum and reached its maximum extent following the Little Ice Age, ∼1,000 years ago (Pelletier et al., 2017). Permafrost conditions strongly influence peatland hydrology and thus the hydrological connectivity of lakes to their surrounding landscape in this region (Quinton et al., 2019). Widespread permafrost precludes deeper groundwater flow but also blocks shallow flow paths and much peatland runoff enters hydrologically isolated depressions rather than connecting with lakes (Quinton et al., 2019;Vitt et al., 2000). Consequently, with less permafrost in warmer climates, there is increased connectivity of lakes with the surrounding peatlands, while lake position in the landscape determines the potential for regional groundwater recharge or discharge (Hokanson et al., 2019). In this study, we assess the influence of climate and permafrost thaw on CO 2 and CH 4 emissions from peatland lakes throughout boreal western Canada. We estimate the current magnitude of greenhouse gas emissions from lakes in the region and use a spacefor-time approach to evaluate how emissions will change with climate warming.

Sampling Sites
The studied region includes Taiga Plains and Boreal Plains ecozones of the Interior Plains in boreal western Canada, an area that extends from 70°N to 56°N. This region encompasses the lands of the Metis, Dene Tha', Woodland Cree, Big Stone Cree in Alberta and the lands of the Inuvialuit Settlement Region, Gwich'in Settlement Area, and the Dehcho of the Northwest Territories. The region is characterized by thick peat deposits and generally flat topography (Hugelius et al., 2020;Zoltai, 1993). The region was previously glaciated and is dominated by glacial and glaciolacustrine surficial geologies (Geological Society of Canada, 2014). Across the latitudinal transect, the ecosystems include similar tree species including aspen (Populus tremuloides), balsam poplar (Populus balsamifera), white spruce (Picea glauca), jack pine (Pinus banksiana), larch (Larix laricina), and black spruce (Picea mariana).
We sampled 20 lakes located in four sites within the study region (Figure 1), spanning across different climatic and permafrost zones (Table S1 in Supporting Information S1). The southern maximum extent of permafrost in the region is approximately 100 km south of the sporadic permafrost zone sampling location (Brown et al., 2002). Projected losses of permafrost thaw have the southern maximum extent of permafrost moving north of the sporadic and discontinuous zones by 2100 under RCP 4.5 and complete near-surface permafrost loss is projected for the study region by 2100 under RCP 8.5 . At each site, we chose five study lakes that had surface areas between 0.3 and 10.1 has (Table S2 in Supporting Information S1), ensuring similar ranges in lake sizes within each location. We only selected lakes surrounded by peatlands and excluded lakes near dirt roads or rivers, which may influence the chemistry of the lakes. Locations of study lakes within the Interior Plains of western Canada. The Interior Plains study region is outlined with black lines, and is comprised of the Taiga Plains (north) and Boreal Plains (south, largely permafrost free) ecozones. Permafrost zones, indicated by blue shadings, are from Brown et al. (2002). Five lakes were monitored for seasonal greenhouse gas emission at each study location, indicated by orange circles. The distribution of peatlands within both the study region and the broader circumpolar region is shown in the top right inset, showing that the study region is one of the major northern peatland regions (Olefeldt et al., 2021). Bottom right photos show one lake from each study location, with all lakes surrounded by treed peatlands.
Sampling took place at all sites during the ice-free seasons in 2018 and 2019. We designed the sampling campaigns to capture seasonal patterns at each site from just after ice-out to just before ice-on. Each lake was sampled a minimum of 7-12 times to get a statistically accurate representation of annual diffusive greenhouse gas emissions (Wik, Thornton, et al., 2016). At the northernmost site, each lake was sampled 7-9 times due to logistical limitations of getting to the sites. Exact sampling dates can be found in the available metadata and in Table S3 in Supporting Information S1.

Water Chemistry
Surface water chemistry was measured from a boat at the location of flux measurements and at the center of the lakes. At each lake, we sampled dissolved CH 4 concentrations using a headspace method (Karlsson et al., 2013) and measured dissolved oxygen, specific conductivity, pH, and temperature in the top 10 cm using a handheld water monitor (Yellow Springs Instrument ProDO). Water samples for DOC and total dissolved nitrogen (TDN) were sampled from the center of the lake following protocols outlined in Burd et al. (2018). Concentrations of DOC andTDN were measured by a TOC-L combustion analyzer with a TNM-L module (Shimadzu, Kyoto, Japan). Specific ultraviolet absorbance (SUVA) at 254 nm was measured from 200 to 700 nm (UV-1280, UV-VIS Spectrophotometer, Shimadzu Corporation, Japan) and was corrected by Milli-Q water blanks. We calculated SUVA by first correcting decadal absorbance at 254 nm (cm −1 ) for interference by Fe as outlined by Weishaar et al. (2003).
Water samples for ion concentrations were filtered and frozen on the day of sampling and kept frozen until analysis. Anions were measured using Ion chromatography (Dionex Ion chromatography DX 600, Thermo Fisher Scientific, US). Cations were measured by ICP-OES (iCAP6300 Duo, Thermo Fisher Scientific, US). Samples for Chlorophyll-a analysis were collected using opaque Nalgene bottles. Duplicate samples for each lake were refrigerated and vacuum filtered on the day of collection using 47-mm diameter, ethanol-rinsed GF/F filters (Whatman). The filters were then placed in petri dishes, wrapped in aluminum foil, and were frozen immediately until analyses at the University of Alberta Biogeochemical Analytical Service Laboratory.

Microbial Sampling
Between July 28 and August 31, 2019, we collected three surface layer (∼15 cm) sediment samples from the edge of each lake where flux measurements were taken. Samples were pooled per lake and DNA was extracted in duplicate from 0.65 g dry sediment per sample using the DNeasy PowerSoil kit (QIAGEN N.V., Hilden, Germany). Sequencing libraries were prepared from 250 ng DNA per sample with the NEBNext Ultra II FS DNA Library Prep kit according to the manufacturer's instructions (New England Biolabs, Ipswich, MA, USA). Samples were sequenced paired end (300 cycles) using a NextSeq 500/550 High Output Kit v2,5 (Illumina, Inc, San Diego, CA, USA). Contamination was negligible, as assessed with both negative (nuclease-free water) and positive (ZymoBIOMICS Microbial Community Standard, Zymo Research Corporation, Irvine, CA, USA) controls. We estimated the abundance of functional marker genes involved in methanogenesis: methyl coenzyme A reductase (mcrA) and methanotrophy: methane monooxygenase subunit A (mmoX) and particulate methane monooxygenase (pmoA), expressed relative to the total number of SILVA 16S hits for each sample with GraftM (Boyd et al., 2018). Whole-metagenome functional profiling relative to the SEED database (used to identify decaheme cytochrome genes and methanogenesis enzyme composition data) was performed with SUPER-FOCUS (Silva et al., 2016), and reads were normalized using the geometric mean of pairwise ratios (Chen et al., 2018). We were unable to extract enough DNA samples from one lake in the discontinuous zone (lake W1) to perform analysis.

Greenhouse Gas Fluxes
Diffusive CH 4 and CO 2 fluxes were measured with small, lightweight chambers (volume = 13.55 L, area = 0.062 m 2 ), covered in reflective tape. Diffusive measurements were taken on the edge of the lakes 1-2 m from the shore to avoid littoral vegetation. For lakes that had unstable edges to deploy the chambers from, we placed the chambers on the water surface from an inflatable boat as to not disturb the sediments. In the few cases when sediment disturbance was observed (i.e., change in the linearity of methane flux and/or higher than normally observed CH 4 concentrations), measurements were stopped and restarted in nearby unaffected locations. We neither avoided nor targeted lake edges with ongoing thermokarst. Most lake edges in the permafrost zones showed evidence of moderate permafrost thaw, but mass wasting of large peat blocks was very uncommon. While there is evidence for both under-and over-estimating lake fluxes based on near shore measurements (Matveev et al., 2016;Schilder et al., 2016), all of the study lakes were surrounded by black spruce trees near the lake edge, providing similarly sheltered wind conditions and comparable environments for comparison across the transect. The chambers were connected to a portable greenhouse gas analyzer (Los Gatos Research, California) in a closed path that recorded the linear change in GHG concentrations within the chamber over time. On each sampling occasion, 3-4 chamber flux measurements of 1-5 min durations (average = 3 min) were taken and then averaged to estimate daytime flux estimate. Typically, each lake was measured on 1-3 days per trip and 7-12 days total. Fluxes were calculated using linear regression (slope) of the change in GHG gas concentration (ppm) inside the chamber over time and the ideal gas law, air temperature, and atmospheric pressure at the time of sampling (Equation 1). Negative slopes for CO 2 measurements were interpreted as CO 2 uptake. CH 4 slopes with R 2 values less than 0.90 were excluded in the daytime average slope calculations.
where E S is the slope from the linear change in gas concentrations within the chamber over time (atm min- T is the average air temperature inside the chamber (K), A is the area of the chamber (m −2 ), M is the molar mass of the CH 4 or C-CO 2 (mol mg −1 ) and E C is a conversion factor to convert from per minute to per day.
We measured ebullitive CH 4 and CO 2 fluxes using a bubble trap (area = 0.073 m 2 ) equipped with a syringe and stopcock. Three out of five lakes at each site were equipped with bubble traps, including three traps along the lake edge and three traps in the lake center. Water column depths below the traps ranged from 0.5 m (L5/L1 edges) to 2.5 m (W1/F2 centers; Table S2 in Supporting Information S1). Traps were set at the beginning of each trip and were sampled 2-4 days later. Accumulated gas bubbles were collected manually using syringes. To measure GHG concentrations in the bubbles we reset the traps and disturbed the sediments, immediately collecting 30 ml of sample, stored in 20 ml glass vials sealed with rubber septa (see in Supporting Information S1 for more discussion on ebullition measurements). Samples were analyzed within 2 weeks after collection at the University of Alberta using an SRI-8610C gas chromatograph equipped with an FID detector and methanizer to convert CO 2 into CH 4 . We calculated the ebullitive flux as outlined by Wik et al. (2013).

Upscaling
We estimated current lake CH 4 and CO 2 emissions from the Taiga Plains by multiplying the area of small and midsized peatland lakes from the Boreal Arctic Wetland and Lake Databset ("BAWLD"; Olefeldt et al., 2021) by the mean CH 4 and CO 2 flux specific to the MAAT reported for each 0.5 × 0.5° grid cell of the database and the number of ice-free days. Daytime diffusive CH 4 emissions were multiplied by a diel correction factor (0.70) to account for diel patterns in CH 4 emissions between day and night (i.e., "24-hr fluxes"; Sieczko et al., 2020). Mean CH 4 and CO 2 emissions and the number of ice-free days were calculated based on linear relationships between mean annual temperature, CH 4 or CO 2 flux, and ice-free days observed across the latitudinal transect ( Figure S4 in Supporting Information S1). Total mean CH 4 emission (corrected diffusion + ebullition) was used in upscaling. We used error propagation of the flux uncertainty and ice-cover day uncertainty to calculate the low and high ranges of current and future total annual carbon emissions. Annual carbon emission values do not include lake area uncertainty from BAWLD, which is estimated to be between −21% and +31% for the low-end and high-end confidence bounds (see Section 4 for more discussion on sources of uncertainty). Carbon dioxide emissions were not corrected for nighttime-factor, but evidence suggests diel patterns in CO 2 respiration are less, bordering negligible for darker, DOC-rich waters (Gómez-Gener et al., 2021), so diel patterns are likely minimal for the lakes that are positive CO 2 emitters. Furthermore, daylight hours across the transect are long, with closer to full days of sunlight throughout most of the summer period in the northern sites. Therefore, applying the diel factor for daily CH 4 emissions is likely a conservative estimate of daily net CH 4 emissions. Though we do acknowledge that the lakes that experience net uptake of CO 2 during the daytime may be small sinks at nighttime-especially for periods closer to fall and spring. We only scaled across grid cells within the Taiga Plains that were within the temperature MAAT range of the study and not the entire extent of the Interior Plains. See Data Availability for information on the land cover database.

Radiative Forcing Model
We used a radiative forcing model that applies simple impulse-response functions and assumed the same warming potential timeframes for CH 4 and CO 2 (Joos et al., 2013). We modeled the carbon pools for CH 4 and CO 2 and mixed atmosphere assumptions following the methods and code by Günther et al. (2020).
While CH 4 was modeled as one pool with a simple exponential decay rate, CO 2 was split into five pools with varying decay rates (Frolking et al., 2006) and did not consider the climatic effects of CO 2 from CH 4 oxidation (Myhre et al., 2013). We compared the radiative forcing trajectories under two global temperature change scenarios (RCP 2.6 and RCP 4.5; IPCC's AR5). To account for the effects of Arctic temperature amplification observed across the Northwest Territories, we multiplied the projected changes in average global temperatures for each RCP by three, which is more conservative than the most recent analysis by the Government of the Northwest Territories which suggest that the region could be warming four times faster (2018). We did not include RCP scenario 8.5 because the projected temperature increases are greater than the range of temperatures experienced across the study transect (range = 9°C; Table S1 in Supporting Information S1). Prior to 2020, we performed a simplified model spin-up, starting in 1900. In the spin-up, we applied current CO 2 and CH 4 fluxes at a constant rate over a 120-year period ( Figure S6 in Supporting Information S1). The purpose of the simplified model spin-up is to bring the model close to an equilibrium point (Turetsky et al., 2020). We then used the year 2020 as a reference point for changes in instantaneous net radiative forcing (Helbig et al., 2017). The model accounts for future changes in CO 2 and CH 4 emissions and the number of ice-free days associated with warming ( Figure S5 in Supporting Information S1). The effect of hotspot emissions along thermokarst lake edges was not explicitly accounted for but was implicitly included since our lake edge sampling included moderately affected thermokarst lake edges. Any influence of permafrost thaw on spring ice-out emissions could not be accounted for and is not considered in the model. Uncertainty for the radiative forcing model was calculated based on one standard deviation of the scaling model coefficients, including flux and ice-free day uncertainty.

Statistical Analysis
All statistical analyses were performed in R statistical software (Version 3.8, R Foundation for Statistical Computing, www.r-project.org). Predictor variables for all models were log-transformed for statistical analyses except water temperature, DOC, and pH. Daytime GHG emission comparisons between zones were performed using mixed-effects models (R Package 3.3.3; Lme4 Package; Bates et al., 2020) with the fixed effect of Permafrost Zone and random effects of the individual lake to account for repeated measures, and season to account for variation in emissions between early, mid, and late-season emissions (Table S4 in Supporting Information S1). Each mixed model for CO 2 , diffusive CH 4 , and ebullitive CH 4 was compared to the respective null model, which only included the random effects, using ANOVA. Individual pairs were then compared using a post-hoc, non-parametric t-test (Bonferroni adjustment; Table S5 in Supporting Information S1). Relationships between environmental variables and GHG emissions were first assessed using a Pearson's correlation matrix ("corrplot" package; Figure S1 in Supporting Information S1). For simple linear regression, structural equation modeling (SEM), and redundancy (RDA) analysis, we used mean values over the entire sampling period and excluded three outlier lakes from analysis in the SEM. The outlier points include lakes L2, L5, and W4 (see Table S2 in Supporting Information S1 for more information on lake IDs). The lakes were designated as outliers due to extreme differences in nutrient concentrations, solute chemistry, and CO 2 emissions, respectively ( Figure S7 in Supporting Information S1). We did not include sulfate in the RDA or further analyses based on the correlation plot which indicated that sulfate covaried with the other groundwater proxies, including base cations and conductivity, and did not have a relationship with CH 4 ( Figure S1 in Supporting Information S1). The RDA was performed using the vegan package's rda() function with CO 2 and diffusive CH 4 emission set as the response variables. Methane and CO 2 fluxes were log transformed to normalize the data. Finally, changes to the functional composition of microbial communities were tested using permutational multivariate analysis of variance with vegan's adonis2 function based on a Euclidean distance matrix.
We used SEMs to quantify the direct and indirect effects of MAAT on environmental variables and average daytime CO 2 and CH 4 emissions (R Package, piecewiseSEM; code adapted from St Pierre et al. (2019). The SEM defines causal links between variables of interest through a system of predefined interconnected linear equations (Shipley, 2009). It then tests the significance of individual paths and the entire model as a whole. The significance of individual paths is determined by a P value < 0.05. However, the model as a whole is defined as a good fit if P > 0.05 from a Fisher's test of directed separation, as this means that the proposed model structure adequately captures the variability in the model variables and the designated paths (Shipley, 2009). Designated model paths must be supported by known links between variables based on prior work. For CO 2, the links and associated references are (a) Chla (as a proxy for productivity) and CO 2 (Engel et al., 2020), (b) phosphate and lake productivity (Carpenter et al., 2001), (c) DOC and productivity (Kissman et al., 2013), and (d) MAAT and phosphate and DOC (Laudon et al., 2012;Roehm et al., 2009;Weyhenmeyer et al., 2015). For CH 4 , the links and prior work included (a) water temperature and CH 4 flux (DelSontro et al., 2016; Yvon-Durocher et al., 2014); (b) water temperature and lake depth (Fang & Stefan, 1999); and (c) MAAT and water temperature (Schneider & Hook, 2010). Iron was not included in the CH 4 SEM because clear links between higher dissolved iron in the continuous permafrost zone could not be established in the prior literature and theories. A summary of the statistical tests used in relation to GHG emissions can be found in Table S6 in Supporting Information S1.

Results and Discussion
We found that lakes in regions with less permafrost and warmer climates had lower CO 2 emissions, or even CO 2 uptake, and much higher CH 4 emissions, than lakes further north (Figure 2). Lakes in warmer climates also had higher pH, concentrations of base cations, DOC, phosphate (PO 4 3− ), and Chla. Between-lake Figure 2. Trends in magnitude and variability of lake greenhouse gas emissions across permafrost zones. Emissions of (a) Diffusive CO 2 and (b) CH 4 were measured during the early, mid, and late ice-free season from five lakes each within the continuous (Cont.), discontinuous (Disc.), and sporadic (Spor.) permafrost zones, as well as south of the permafrost region (None). Boxplots indicate the 25th, median, and 75th percentiles of emissions, and the yellow diamonds indicate the average emissions. Different letters above the boxplots for CH 4 indicate statistically significant differences between permafrost zones, accounting for repeated measurements throughout the season (see Section 2 and Table S4 in Supporting Information S1). Diffusive and ebullitive CH 4 flux significant differences are presented separately. Positivie CO 2 fluxes represent a net source to the atmosphere while negative values represent a net sink. Outlier CO 2 fluxes below −800 and above +800 mg C m −2 d −1 are not shown. Methane fluxes have not been corrected for diel variation in this figure.
variability also increased from colder to warmer climates for both CO 2 and CH 4 emissions, and water chemistry ( Figure 2). While CO 2 and CH 4 emissions both varied with latitude, they were only weakly related to each other ( Figure S2 in Supporting Information S1). Below, we explore the distinct drivers of CO 2 and CH 4 emissions and use the information to project future emissions and the net radiative forcing resulting as feedback to climate warming under different RCP scenarios.

Connectivity, Primary Productivity, and CO 2 Emissions
Our study shows that climate change is likely to reduce peatland lake CO 2 emissions over the ice-free period, as its positive effects on primary productivity appeared greater than effects on respiration across the transect. We attribute this result primarily to indirect effects of climate through increased hydrological connectivity of lakes to the surrounding landscape. Across the transect, diffusive CO 2 emissions were lower for lakes with higher concentrations of Chla, which also were lakes with higher DOC, light absorbance, PO 4 3− , and pH ( Figure 3a). These water chemistry attributes are characteristic of near-surface water in nutrient-rich groundwater-connected peatlands (i.e., fens) that are relatively more widespread in discontinuous and non-permafrost zones (Vitt et al., 2000;Wood et al., 2016; Figure 4). While there was a trend of decreasing CO 2 emissions in warmer climates, the SEM highlighted that this link was mediated through effects on lake water chemistry (R 2 = 0.78, Fisher's C = 12.9, P = 0.40; P > 0.05 = good fit; Figure 3b). In the SEM, we hypothesized that higher MAAT allows for less permafrost and increased lake connectivity with groundwater-fed peatlands, raising lake DOC, phosphate, and pH, which in turn increases primary productivity, indicated by higher Chla (Engel et al., 2020), and shifts the carbonate equilibrium from CO 2 to bicarbonate, both causing reduced CO 2 emissions. Notably, while we did not design the study to monitor CO 2 emissions from abrupt thaw features, our study does include lake edges impacted by moderate levels Figure 3. Influence of lake water chemistry and lake characteristics on greenhouse gas emissions across permafrost zones. (a) Redundancy analysis (RDA) ordination for CO 2 and CH 4 emissions, where explanatory variables include water chemistry variables along with maximum lake depth and mean annual air temperature (MAAT). The circles represent individual lakes, colored by permafrost zonation. Close spatial proximity of the circles to one another indicates that lakes are more similar in their emissions, whereas circles that are far away indicate dissimilarities. The direction of the arrows represents the relationship between the environmental variables and CO 2 and CH 4 flux. Variables with opposite facing arrows are negatively correlated and those with arrows leading in the same direction and orientation are positively correlated. Orthogonally related arrows indicate no relationship. The percent of variability explained by each ordination axis is shown in parentheses. The RDA shows that lake CO 2 and CH 4 fluxes were not correlated, with CH 4 emissions associated with lake depth, water temperature, and MAAT while CO 2 emissions were more closely associated with higher concentrations of nutrients, DOC, and chlorophyll a. Base_ Cations represents summed Ca, Mg, and Na concentrations. Potassium (K) is shown separately as it followed a dissimilar trend to the other cations. Structural equation models were constructed individually for (b) CO 2 and (c) CH 4 . Arrows with solid lines and asterisks represent statistically significant pathways (p values = * <0.05, ** <0.01, *** <0.001), while dotted lines represent non-significant pathways. The values by each arrow represent standardized model coefficients, indicating positive and negative relationships. Goodness of fit (Fisher's C values) were 12.90 and 4.85, respectively for the CO 2 and CH 4 models, while P values were 0.40 and 0.30. P values > 0.05 indicate that model structure does not differ from that expected by the data. Outlier lakes were excluded from the RDA and SEM analysis (see Section 2.6).
of thaw. Thus we may have missed small hotspots of higher CO 2 emissions from time-sensitive abrupt thaw locations, however, these features were rare across the transect and we consider our sampling locations to be representative of lake shorelines in the region. CO 2 emission from ebullitive fluxes did not follow a strong north to south trend but emissions were generally higher from the non-permafrost zone than the continuous zone ( Figure S3 in Supporting Information S1), suggesting a potential link between temperature and sediment respiration. However, the magnitude of ebullitive CO 2 fluxes was minor compared to diffusive emissions (∼7% of the total flux), and did not change the latitudinal trends or net sink/source capacity of individual lakes.
Our findings of decreasing diffusive CO 2 emissions from lakes with higher DOC concentrations contrast most other boreal lake surveys (Lapierre et al., 2015;Weyhenmeyer et al., 2015). Previous studies have shown strong positive relationships between lake DOC concentrations and CO 2 emissions, suggesting that greater connectivity to terrestrial sources of DOC either enhances lake respiration and/or suppresses primary productivity through light attenuation, resulting in higher net CO 2 emissions (Ask et al., 2012;Lapierre et al., 2015). We speculate that the difference between our results and previous studies stems from our inclusion solely of small peatland lakes with relatively high DOC concentrations (range 17-50 mg L −1 ). Other surveys include lakes of various origins and sizes, and thus lakes with much lower DOC concentrations (Ask et al., 2012;Lapierre et al., 2015). For the peatland lakes in this study, it is possible that further elevated DOC concentrations in warmer climates have limited additional influence on respiration, though notably, we did not directly measure respiration or metabolism. However, our results do suggest that higher DOC may promote primary productivity through the delivery of complex nutrients. For example, while the SEM showed a non-significant direct link between DOC and CO 2 flux, there was a significant causal link between DOC and Chla, indicating DOC influences primary productivity and thus indirectly impacts net CO 2 Figure 4. Influence of climate and permafrost conditions on greenhouse gas emissions from peatland lakes. Lakes in the continuous permafrost zone were found to be net sources of CO 2 and had low CH 4 emissions, while lakes outside the permafrost region were CO 2 sinks and had high CH 4 emissions. Absence of permafrost allows for increased hydrological connectivity to the surrounding landscape, potentially both deeper flow-paths through mineral soils and near-surface flow-paths through the surrounding peatlands. In this study, we found that lakes outside the permafrost region had generally higher pH and higher concentrations of ions, nutrients, and chlorophyll-a-suggesting that elevated primary productivity (PP) outweighs respiration (R) in lakes outside the permafrost region. Lake diffusive CH 4 emissions were primarily linked to lake water temperatures and were 8-times greater south of the permafrost boundary compared to within the continuous permafrost zone. emissions (Figure 3b). Our results suggest that peatland lake CO 2 emissions are likely to respond to climate change in a different way than other lake types with varying DOC and nutrient conditions, such as deeper lakes in regions of Canada and Scandinavia underlain by Precambrian rock, and may also vary regionally among peatland lakes depending on landscape characteristics.

Temperature Sensitivity, Microbial Dynamics, and CH 4 Emissions
Both ebullitive and diffusive CH 4 emissions were greater in warmer climates, with emissions from lakes in the non-permafrost zone being ∼8 times greater than emissions from lakes in the continuous permafrost zone throughout the ice-free season. Between-lake variability in CH 4 emissions was also greatest in the sporadic and non-permafrost zones (Figure 2b). Ebullitive emissions represented ∼28% of the total C emissions from lakes where it was measured; a fraction that did not vary systematically across the transect. Ebullitive CH 4 emissions were generally highest in spring, while diffusive emissions were highest during summer or fall ( Figure S3; Table S5 in Supporting Information S1). Higher MAAT and water temperature were strongly related to higher diffusive CH 4 emissions, while higher base cation concentrations, decreasing lake depth, and lower concentrations of iron were weakly associated with higher emissions (Figure 3a). The lack of a strong relationship between diffusive CH 4 emissions and lake size and depth, as is often found (Holgerson & Raymond, 2016;Wik, Varner, et al., 2016), was not surprising as all studied lakes were relatively small, shallow, and predominately non-stratified (Table S2 in Supporting Information S1). The SEM analysis showed a significant relationship of MAAT with water temperature, which in turn influenced diffusive CH 4 emissions (R 2 = 0.67, Fischer C = 4.3, P = 0.3; Figure 3c). The apparent temperature sensitivity (Q 10 ) of diffusive CH 4 emissions across the transect to changes in water temperature was 16, substantially higher than the Q 10 of ∼4 that is commonly found for the direct temperature response of CH 4 emissions within individual ecosystems (Yvon-Durocher et al., 2014).
The high apparent temperature sensitivity of CH 4 emissions across the transect suggests the observed direct response to water temperature was amplified by additional factors, possible shifts in microbial communities, redox conditions, or substrate quality (Rasilo et al., 2015). Of these, we found the least support for increasing CH 4 emissions being driven by shifts in substrate quality. Greater lake primary productivity in warmer climates was hypothesized to increase the availability of labile substrate for methanogenesis (Sepulveda-Jauregui et al., 2018), but we found only a weak relationship between CH 4 emissions and Chla ( Figure 3a). However, we collected Chla from the lake centers which may not capture primary productivity related to littoral vegetation and associated root exudates. Permafrost thaw was also hypothesized to increase the availability of labile substrates through thermokarst lake expansion. However, visual evidence of ongoing thermokarst lake expansion was found in all permafrost zones, and the highest CH 4 emissions were found in the non-permafrost zone. Lastly, it is possible that shifts in surrounding wetland and littoral vegetation toward denser, taller emergent macrophyte species in warmer climates could provide external sources of labile organic substrates and dissolved CH 4 to the lake (Emilson et al., 2018), although we did not observe any apparent links between shoreline vegetation and CH 4 emissions.
Several lines of evidence instead suggested that shifts in redox conditions and microbial communities contributed to the high climate sensitivity of CH 4 emissions. Our assessment of microbial functional gene abundances showed that warmer climates were associated with a preferential increase of the methanogenesis marker methyl coenzyme A reductase (mcrA) over the methanotrophy markers, methane monooxygenase subunit A (mmoX), and particulate methane monooxygenase (pmoA) ( Figure 5). These findings correlate with observations of greater increases in CH 4 production compared to CH 4 consumption in lake sediments with high production rates (D'Ambrosio & Harrison, 2021). While gene abundances do not necessarily correlate with microbial activity (Freitag & Prosser, 2009), this trend was similar to that of observed CH 4 emissions (Figure 2b). A lower abundance of mcrA in the northern lakes may be reflective of a temperature effect, wherein methanogenesis is suppressed under colder conditions. However, we also found evidence suggesting suppression of methanogenesis in the continuous permafrost zone through greater dominance of more thermodynamically favorable iron reduction. Lakes in the continuous permafrost zone trended towards higher concentrations of iron (Fe) in the water (Table S2 in Supporting Information S1), and higher abundances of Fe(III) reduction marker genes (decaheme cytochromes) in sediments ( Figure  S4 in Supporting Information S1). These conditions could also promote anaerobic CH 4 oxidation (Ettwig et al., 2016). Higher normalized counts of genes involved in iron acquisition and metabolism in colder climates were furthermore associated with a shift in enzymatic composition from hydrogenotrophic and acetoclastic methanogenesis towards the stoichiometrically less productive methylotrophic form (Rhew et al., 2007; Figure S4 in Supporting Information S1). The mechanisms driving higher dissolved iron and abundances of decaheme cytochromes within the continuous permafrost zone are unclear, but the negative relationship between lake iron, pH, and base cation concentrations suggests that hydrological connectivity to deeper groundwater sources is important (Figure 3a).
Our study suggests that peatland lake CH 4 emissions have a high climate sensitivity, likely driven by increasing CH 4 production in organic-rich sediments. High climate sensitivity of peatland lakes is noteworthy as their current CH 4 emissions are higher than emissions from most other boreal lake types . We note however that the observed CH 4 emissions in this study (average diffusion = 27 mg CH 4 m −2 d −1 and ebullition = 11 mg CH 4 m −2 d −1 ) were in the low range of what has been reported previously for boreal peatland lakes (diffusive interquartile range [IQR] = 9-101 mg CH 4 m −2 d −1 , ebullitive IQR = 3-89 mg CH 4 m −2 d −1 ; Kuhn et al., 2021). Furthermore, our latitudinal trends in GHG emissions were more pronounced than those found from lakes in other northern regions including other peatland regions (Sepulveda-Jauregui et al., 2015;Serikova et al., 2019). Differences in lake CH 4 emissions between major peatland regions may be due to landscape characteristics and Holocene history which influence water chemistry and organic matter quality, but differences can also be due to climate and measurement methods (including timing, placement, and number of measurements). Our study implements a methodology that constraints ice-free season, diffusive, and ebullitive CH 4 fluxes from a large number of lakes across a permafrost transect-yielding robust findings.

Future Peatland Lake Emissions and Net Radiative Forcing
We used an atmospheric perturbation radiative forcing model to assess the relative importance of projected shifts in peatland lake CO 2 and CH 4 emissions within the Taiga Plains region until 2100 (Frolking et al., 2006). This model considers the larger radiative forcing efficiency of CH 4 compared to CO 2 (Myhre et al., 2013), and the shorter atmospheric lifetime of CH 4 (Allen et al., 2018). We assumed that emissions and the number of ice-free days will change linearly with increasing temperatures, based on the trends observed across our transect ( Figure S5 in Supporting Information S1), and assessed changes in emissions based on the RCP Scenarios 2.6 and 4.5, adjusted for temperature amplification in the Taiga Plains (Figure 6a). This model does not include potential emissions from abrupt permafrost thaw and drastic land slumping, which is projected to have an intense, but short-lived, impact on net GHG emissions from lakes (Walter Anthony Figure 5. Trends in methane-related functional genes among permafrost zones. (a) Relative abundance of the (a) Methanogenesis marker mcrA (methyl coenzyme A reductase) and methanotrophy markers (b) pmoA (particulate methane monooxygenase) increased and (c) mmoX (methane monooxygenase subunit-A) moving south along the transect. Boxplots represent median and quartile ranges for samples taken from each lake once mid-summer 2019 from the continuous (Cont.), discontinuous (Disc.), and sporadic (Spor.) permafrost zones, as well as south of the permafrost region (None). Note the different y-axis scales. Differences are not statistically significant, likely due to low sample sizes (n = 5 for every zone except the discontinuous zone where n = 4). et al., 2018). Thus, our results may under-estimate net radiative forcing for both CH 4 and CO 2 , however, our study design does take into account emissions from moderately impacted thermokarst edges which were the most common thermokarst feature across the study lakes. Peatland lakes <10 km 2 in size cover ∼260,000 km 2 within the Taiga Plains (Olefeldt et al., 2021) and were estimated to emit 0.20 ± 0.08 Tg CO 2 -C y −1 and 0.048 ± 0.01 Tg CH 4 -C y −1 under the current climate. Emissions of CO 2 were estimated to decrease, on average, by 16.5% and 68% under the RCP 2.6 and 4.5 scenarios, respectively (Figure 6b), while CH 4 emissions increased by 31% and 121%, respectively (Figure 6b). Overall, increasing CH 4 emissions dominated the net radiative forcing under both RCP scenarios. Increasing CH 4 emissions added, on average, 5 and 15 fW m −2 by 2100 under RCP 2.6 and 4.5, respectively, which was only partially offset by reduced CO 2 emissions equivalent to −1 and −2 fW m −2 , respectively, when compared to the no-warming scenario (Figure 6c).

Conclusion
Our study shows that climate change is likely to cause reduced CO 2 emissions but greatly increased CH 4 emissions from boreal peatland lakes in western Canada over ice-free periods (Figure 3). Our projections of future peatland lake CO 2 and CH 4 emissions are limited to the ice-free season, and rely on a space-for-time approach with inherent assumptions, but are likely to be robust. Our annual GHG emission estimates contain many uncertainties including measured fluxes and variability in ice-free days across the transect. Further data are required to assess ice-out emissions in the region and associated uncertainties. Several studies have shown that ice-out emissions generally scale with emissions from the ice-free season (Wik et al., 2014), however, ice-out emission studies across latitudinal gradients and within the context of potential changes Figure 6. Net radiative forcing due to changing CO 2 and CH 4 emissions from Taiga Plains peatland lakes under RCP 2.6 and 4.5 scenarios. (a) The global RCP 2.6 and 4.5 temperature scenario pathways for the 21st century have been multiplied by a factor of three to account for Arctic amplification in the Taiga Plains, western Canada (Government of the Northwest Territories, 2018). (b) Modeled current and future CO 2 and CH 4 emissions and total emissions (CO 2 + CH 4 ) from peatland lakes on the Taiga Plains. Error bars represent 95% confidence intervals from the model output which includes uncertainty from fluxes and the number of ice-free days. (c) The 21st century radiative forcing resulting from altered lake CO 2 and CH 4 fluxes under the RCP 2.6 and 4.5 scenarios. Radiative forcing is referenced to the year 2020. A scenario of no warming is included, under which the net radiative forcing is only influenced by lake CO 2 emissions while steady CH 4 emissions are in equilibrium with the atmosphere. Radiative forcing due to altered lake CO 2 and CH 4 emissions is expressed as 10 −15 fW divided by the total lake area of the Taiga Plains. The shaded areas represent the standard deviations from the flux models. Under both the RCP 2.6 and 4.5, we find that the net radiative forcing is dominated by increasing CH 4 emissions and only partly offset by reduced CO 2 emissions. in net annual GHG emissions with climate change are still lacking. Additionally, most of the existing iceout emissions are based on comparing water column and CH 4 storage before and after ice-out, and the knowledge about how much of the stored CH 4 is emitted versus oxidized is still limited. Our annual and future GHG estimates do not include uncertainty around the lake area. Additional uncertainty associated with small and midsize peatland lakes in BAWLD is estimated to be between −21% and +31%. We expect the true annual emissions to vary within this range, but it is most likely that small lakes are under-counted (Messager et al., 2016;Olefdt et al., 2021), implying our annual and future estimates are likely conservative. Furthermore, we acknowledge that our measurements represent daytime conditions and that nighttime conditions may lead to higher respiration in some lakes, and may alter the net ice-free estimates of CO 2 exchange. However, our results still highlight key mechanistic differences between CO 2 and CH 4 dynamics in lakes across the study region. Water chemistry and microbial community data allowed for the identification of processes likely responsible for the observed, and opposing CH 4 and CO 2 emissions trends along the transect. We showed that peatland lake CO 2 and CH 4 emissions were driven by different processes and thus not strongly inter-correlated. Both the direct effect of warming and the indirect effect of permafrost thaw was important, where the absence of permafrost in warmer regions allowed for greater variability in lake hydrological connectivity and thus also greater variability in CO 2 and CH 4 emissions. The influence of permafrost thaw on lake productivity and reduced CO 2 emissions, and the extremely high-temperature sensitivity of CH 4 emissions, are likely associated with characteristics of peatland lakes and the hydrology of their surrounding landscape. While additional data from other major peatland regions are needed to support our findings, our study strongly suggests that peatland lakes, which make up ∼18% of all northern lake areas (Olefeldt et al., 2021), need to be considered separately from other lake types when assessing the impact of climate change on future GHG emissions.