Climate warming has direct and indirect effects on microbes associated with carbon cycling in northern lakes

Northern lakes disproportionately influence the global carbon cycle, and may do so more in the future depending on how their microbial communities respond to climate warming. Microbial communities can change because of the direct effects of climate warming on their metabolism and the indirect effects of climate warming on groundwater connectivity from thawing of surrounding permafrost, especially at lower landscape positions. Here we used shotgun metagenomics to compare the taxonomic and functional gene composition of sediment microbes in 19 peatland lakes across a 1600‐km permafrost transect in boreal western Canada. We found microbes responded differently to the loss of regional permafrost cover than to increases in local groundwater connectivity. These results suggest that both the direct and indirect effects of climate warming, which were respectively associated with loss of permafrost and subsequent changes in groundwater connectivity interact to change microbial composition and function. Archaeal methanogens and genes involved in all major methanogenesis pathways were more abundant in warmer regions with less permafrost, but higher groundwater connectivity partly offset these effects. Bacterial community composition and methanotrophy genes did not vary with regional permafrost cover, and the latter changed similarly to methanogenesis with groundwater connectivity. Finally, we found an increase in sugar utilization genes in regions with less permafrost, which may further fuel methanogenesis. These results provide the microbial mechanism for observed increases in methane emissions associated with loss of permafrost cover in this region and suggest that future emissions will primarily be controlled by archaeal methanogens over methanotrophic bacteria as northern lakes warm. Our study more generally suggests that future predictions of aquatic carbon cycling will be improved by considering how climate warming exerts both direct effects associated with regional‐scale permafrost thaw and indirect effects associated with local hydrology.

climate warming on their metabolism and the indirect effects of climate warming on groundwater connectivity from thawing of surrounding permafrost, especially at lower landscape positions. Here we used shotgun metagenomics to compare the taxonomic and functional gene composition of sediment microbes in 19 peatland lakes across a 1600-km permafrost transect in boreal western Canada. We found microbes responded differently to the loss of regional permafrost cover than to increases in local groundwater connectivity. These results suggest that both the direct and indirect effects of climate warming, which were respectively associated with loss of permafrost and subsequent changes in groundwater connectivity interact to change microbial composition and function. Archaeal methanogens and genes involved in all major methanogenesis pathways were more abundant in warmer regions with less permafrost, but higher groundwater connectivity partly offset these effects. Bacterial community composition and methanotrophy genes did not vary with regional permafrost cover, and the latter changed similarly to methanogenesis with groundwater connectivity. Finally, we found an increase in sugar utilization genes in regions with less permafrost, which may further fuel methanogenesis. These results provide the microbial mechanism for observed increases in methane emissions associated with loss of permafrost cover in this region and suggest that future emissions will primarily be controlled by archaeal methanogens over methanotrophic bacteria as northern lakes warm. Our study more generally suggests that future predictions of aquatic carbon cycling will be improved by considering how climate warming exerts both direct effects associated with regional-scale permafrost thaw and indirect effects associated with local hydrology.

K E Y W O R D S
aquatic carbon cycling, community composition, functional gene analysis, groundwater connectivity, metagenomics, methanogenesis, methanotrophy, permafrost thaw

| INTRODUC TI ON
Northern lakes are disproportionately important in the global carbon cycle (Hastie et al., 2018;Wik et al., 2016), but there is uncertainty about how they will respond to a warming climate. Much of the carbon that is processed in lakes occurs in sediments where organic matter accumulates and is processed by microorganisms (Bastviken et al., 2004;Gudasz et al., 2010). Changes in microbial community composition can consequently dictate how lake carbon cycling responds to climate warming (Glassman et al., 2018;Schimel & Schaeffer, 2012). Microbial processing of carbon is also highly temperature dependent (Yvon-Durocher et al., 2014;Zhu et al., 2020), and some processes may be more sensitive to climate warming than others . For example, laboratory incubations suggest that as mineralization of sediment organic matter by methanogenic microbes increases, consumption of methane by methanotrophic organisms may become increasingly inefficient (D'Ambrosio & Harrison, 2021), but supporting field data are lacking. Therefore, as lake-rich northern latitudes warm faster than the rest of the planet (Rantanen et al., 2022), there may be large but unknown consequences for carbon cycling arising from changes in microbial composition and function.
In addition to directly accelerating metabolic rates, warmer temperatures indirectly influence microbes involved in the carbon cycling of northern lakes by changing local environments. Warming can specifically increase the concentration and diversity of organic substrates available for heterotrophic respiration and change the energetic favourability of reduction-oxidation (redox) reactions (Conant et al., 2011;Johnston et al., 2019). These impacts may be large in northern lakes surrounded by perennially frozen ground known as permafrost, which stores about twice as much carbon as in the atmosphere (Schuur et al., 2015). As permafrost thaws, hydrological connectivity increases (Connon et al., 2014), and lateral flow paths transport previously frozen organic matter downstream (O'Donnell et al., 2012;Wright et al., 2009). Terrestrial organic matter inputs can explain much of the variation in the microbial composition of northern lake sediments (Orland et al., 2019). Once flow paths deepen into the mineral layer they can become enriched in dissolved ions from weathering, including terminal electron acceptors (TEAs) like ferric iron and sulphate, which can further shift the activity of microbes in downstream waters (Lipson et al., 2013;Walvoord et al., 2012). Inorganic TEAs thermodynamically favour respiration pathways other than methanogenesis (Rissanen et al., 2017), and can limit decomposition altogether by adsorbing to organic matter and making it less physically available to microbes (Lalonde et al., 2012). Groundwater from these deepened flow paths can also be discharged through highly organic fens and riparian soils, thereby carrying high concentrations of both TEAs and dissolved organic carbon (DOC) (Frey & McClelland, 2009). Lakes in lower landscape positions will be particularly suited to intercept groundwater discharge after permafrost thaw (Hokanson et al., 2019), and thus potentially have higher concentrations of TEAs and higher pH. In landscapes at higher positions or with more intact permafrost, lakes will lack this groundwater hydrological connectivity (Frey & McClelland, 2009). The complete loss of permafrost, especially in regions with excess ground ice, can contrast the inhibitory effects of hydrological connectivity on some microbial processes by causing ground subsidence and landscape flooding. The resulting shallower water tables and highly reduced conditions will favour certain reactions, in particular methanogenesis (Hodgkins et al., 2014) and anaerobic methane oxidation (Perryman et al., 2020), and can restructure microbial community composition (Ruuskanen et al., 2018). For these reasons, hydrological connectivity may act as a general indicator of the local microbial processes in permafrost landscapes alongside regional climate gradients.
The functioning of microbial communities may also respond to climate warming independently of taxonomic composition, including in environments affected by permafrost (Frank-Fahle et al., 2014;Liu et al., 2021;Mackelprang et al., 2011). Microbial diversity studies have consistently shown that phenotype, not taxonomy, determines the survival of organisms in a given environment, so they cannot reliably be inferred from one another (Inkpen et al., 2017). For example, Chen, Liu, et al. (2021) found less taxonomic diversity after permafrost thaw, but more functional gene diversity and greater abundances of genes involved in the degradation of xylan and aromatic compounds . However, many studies infer functional changes from 16S metabarcoding Seward et al., 2020), and so cannot discern underlying functional changes in the absence of taxonomic ones. Methanogenesis, for example, is a highly diverse process. Tracking it by amplicon sequencing of the marker gene mcrA (methyl coenzyme A reductase) cannot differentiate among the different metabolic steps involved in methanogenesis and that can vary in abundance across the environment (Dziewit et al., 2015). By contrast, shotgun metagenome sequencing can reveal functional changes that are absent from taxonomic analyses while also providing more robust taxonomic assignments (Clooney et al., 2016;Durazzi et al., 2021). Most field studies employing this approach sample permafrost-affected soils and wetlands (Rivkina et al., 2016;Yuan et al., 2018), but not lake sediments, despite the latter being a key location for carbon cycling (Zandt et al., 2020). Existing permafrost studies also mostly sample local thaw gradients (Hultman et al., 2015;Woodcroft et al., 2018) or the soil profile at a single site (Xue et al., 2020). These approaches miss an opportunity to compare the indirect effects of variation in local hydrology associated with climate warming and direct effects of regional warming on microbial diversity and function (Comte et al., 2016;Ebrahimi & Or, 2017).
Here we tested how the taxonomic and functional composition of microbial communities in lake sediments changes as permafrost cover is reduced. Unlike previous studies that examine thaw gradients in a single catchment (Peura et al., 2020) or after artificially induced thaw (Mackelprang et al., 2011), we sampled 19 lakes across a 1600-km permafrost transect. The transect acted as a space-fortime substitution of future climate warming and allowed us to separate regional versus local impacts of climate warming. All the lakes were surrounded by peatlands within the Interior Plains of western Canada, with similar underlying surficial geology, but were located within different regions of permafrost cover: none, sporadic, discontinuous and continuous. Variation in thaw extent and lake position within each region further allowed us to separate the local, indirect effects of warming mediated through groundwater connectivity from the direct effects of regional climate warming represented by permafrost region. Using shotgun metagenomic sequencing, we related the abundances of both microbial genera and functional genes involved in carbon cycling to environmental conditions. We expected that microbial communities would vary in taxonomic and functional gene composition across the latitudinal gradient with differences in regional climate, for example, Seward et al. (2020). We also expected that groundwater connectivity within sites would change the composition and reduce the total abundance of both methanogenic communities and functional genes involved in organic matter decomposition by increasing the availability of inhibitory non-O 2 TEAs (Miller et al., 2015).

| Study sites
We studied 19 lakes ranging from 0.3 to 10.1 ha in area and distributed between 56.0°N and 67.5°N in boreal western Canada ( Figure 1). Five lakes were sampled in each of three regions of permafrost coverage (continuous, sporadic and none), and four lakes were sampled in the discontinuous permafrost region. The primary difference among regions was in mean annual temperature, which varied between −7.3°C and +1.5°C (Fick & Hijmans, 2017), with no other major bioclimatic gradient. All the study regions were generally dominated by flat topography with thick overburden underlain by erodible sedimentary bedrock that influences groundwater chemistry (Ecosystem Classification Group et al., 2009). All the lake catchments were also dominated by similar treed peatlands; see Kuhn et al. (2021) for further details. The study region encompassed 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.

| Field sampling
We sampled lake sediments and environmental characteristics during late-summer (28 July 2019 to 31 August 2019) when we expected active layers to be thawing. At each lake, we collected the top 15 cm of sediment using a sterile trowel and sterile gloved hands at three locations along the stable non-thermokarst edge of each lake. Samples were immediately pooled in individual sterile bags at a maximum 2-8°C as permitted by field storage facilities, and stored at −20°C within 5 days, and −80°C within 14 days. Concurrently, we measured electrical conductivity (EC) of the overlying water as an indicator of connectivity to mineral groundwater sources and refer to it as such hereafter (Burd et al., 2018;Olefeldt et al., 2013; Wlostowski F I G U R E 1 Location of study sites within western Canada. At each orange circle, we sampled five lakes; four in the discontinuous permafrost region (total n = 19 lakes). All sites were within the Interior Plains physiographic region outlined by black lines. Permafrost regions are indicated in blue from Heginbottom et al. (2002). The dominance of peatlands within the study region and their distribution with the broader circumpolar region is shown in the top right inset . Bottom right photos show one lake from each study region, with all lakes surrounded by treed peatlands. Figure et al., 2016). The measurements were taken from surface water in the centre of each lake using a handheld probe that also recorded pH and water temperature (Yellow Springs Instrument, ProDO).
We also measured EC, pH and water temperature alongside water chemistry on separate days at each site to understand how groundwater connectivity was a general indicator of environmental conditions. Surface waters were collected on 7-12 occasions between May 2018 and October 2019 from the centre of each lake.
Water was passed through glass fibre filters into precombusted glass amber bottles and acidified to pH <2 with HCl to measure DOC concentration on a TOC-L combustion analyser (Shimadzu Corporation).
Ultraviolet (UV) absorbance at 254 nm of filtered but unacidified samples was measured with a UV-1280 spectrophotometer (Shimadzu Corporation), and corrected by Milli-Q water blanks. We divided absorbance by DOC concentration to derive specific UV absorbance (SUVA), which reflects the aromaticity of dissolved organic matter. SUVA was corrected for interference by Fe (Weishaar et al., 2003), which was measured in an additional filtered water sample that was immediately frozen in the field for analysis of ion concen- a DX 600 ion chromatography system (Thermo Fisher Scientific) and iCAP6300 Duo inductively coupled plasma-optical emission spectrometer (Thermo Fisher Scientific). We then averaged all the measurements across sampling dates to derive a long-term average description of each site. Both EC and pH were highly correlated between the late-summer 2019 measurements and long-term means (r = .99 and .92 respectively; p < .001 for both). Sequencing libraries were prepared from 250 ng DNA per sample.

| Sequencing of microbial communities
Library preparation was performed with NEBNext Ultra II FS DNA Library Prep and single barcoding with NEBNext Multiplex Oligos for Illumina according to the manufacturer's instructions (New England Biolabs). Libraries were quantified using the Qubit dsDNA High Sensitivity Assay Kit on the Qubit Fluorometer. Library quality was examined on a Bioanalyzer HS DNA chip (Agilent) before pooling at equimolar concentrations into a single multiplexed sample. Samples were sequenced on a NextSeq 500 using a NextSeq 500/550 High Output Kit v2 (300 cycles, paired-end) (Illumina, Inc.). Composition in the mock community was not statistically different from the theoretical distribution provided by the manufacturer ( Figure S1). The negative control did not yield sufficient DNA to proceed with sequencing, suggesting there were no significant laboratory contaminants. Sequences are available in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB36558 (https:// www.ebi.ac.uk/ena/brows er/view/PRJEB 36558).
Taxonomic assignment of reads (>Q30) to the genus level was performed with kraken2 using the RefSeq PFP Plus October 2020 database (Wood et al., 2019). Functional profiles of samples were grouped into SEED subsystems using SUPER-FOCUS based on the RapSearch2 aligner (Silva et al., 2016;Zhao et al., 2012). SEED subsystems consist of groups of genes that are organized into biological pathways or complexes that undertake different functions (Overbeek et al., 2005). We focused our analyses on functions with the most detail about their biological roles called Subsystem 3.
Taxonomic and functional profiles were then normalized using the geometric mean of pairwise ratios to account for the compositional nature of the data (Chen, Reeve, et al., 2018). We subsequently removed taxa from each sample with a relative abundance of <0.1% based on the accuracy of taxonomic assignment for the mock community positive control ( Figure S1). We repeated our statistical analyses at a lower threshold (0.05%) to confirm that this filtering threshold did not change our conclusions by excluding rare taxa.

| Statistical analyses
We tested if EC reflected groundwater connectivity using redundancy analysis (RDA). We constrained the average seasonal concentrations of DOC and major ions derived from mineral groundwater (Ca, Cl, Fe, Mn, Mg, Na, S, Zn) onto EC of each site using the rda function in vegan. Variables that were not normally distributed were log-transformed and all variables were separately normalized to a mean of zero and standard deviation of one. We performed 999 permutations of the concentration data to test the statistical significance of the RDA.
We partitioned the variance in community composition of all microbes, archaea and bacteria explained by climate, hydrology, organic matter inputs and lake morphometry using RDA with the varpart function from vegan. The climate variables were permafrost region and mean annual temperature and precipitation between 1970 and 2000 from WorldClim v2 (Fick & Hijmans, 2017).
For organic matter inputs, we used DOC concentration and SUVA measured during the longer-term seasonal monitoring in July 2019 and as close as possible to sediment sampling. We similarly used EC measured on these dates to describe groundwater connectivity. By measuring these variables before sediment collections, we could specifically test the responses of microbes to antecedent conditions. EC was also highly correlated between July 2019 and sediment sampling in late July to early August 2019 (r = .98; p < .001), suggesting that conditions persisted throughout the summer. We also expected hydrology would be influenced by the per cent of peatland cover in the surrounding catchment (Drexler et al., 1999). Finally, the lake characteristics included the ratio of lake area to catchment area and maximum lake depth. These variables further allowed us to characterize differences in landscape positions of lakes within each permafrost region. Predictors that were not normally distributed were log-transformed and we used the Hellinger transformation to standardize community composition data. The marginal effect of each partition was further tested using a permutation test that randomized the community data 999 times.
We also tested if more geographically distant sites had more dissimilar microbial communities. We used a Mantel permutation test to correlate a matrix of pairwise geodesic distances between sites with a matrix of their pairwise differences in composition. Compositional dissimilarities were calculated with a Horn-Morisita dissimilarity index (Morisita, 1959) because this index is robust to the effects of sample size and diversity while being appropriate for transformed (normalized) data as used here (Wolda, 1981).
Finally, we tested if the functional genes of microbial communities varied among permafrost regions and with groundwater connectivity. SEED subsystems consist of groups of genes that are organized into biological pathways or complexes that undertake different functions (Overbeek et al., 2005). We focused our analyses on functions with the most detail about their biological roles called Subsystem 3. We focused on SEED Subsystem 3 functions associated with methanogenesis, organic matter decomposition and both sugar utilization and energy production by microbes. We first compared the functional gene composition of each SEED function using a PERMANOVA with Bray-Curtis distances, and visualized functions with statistically significant effects using non-metric multidimensional scaling ordination with three dimensions. We also fitted separate linear models to the log-transformed sum of normalized abundances for all genes in each Subsystem 3 function, with permafrost region and log-transformed EC as predictors. For these analyses, we generated a fourth category of functions associated with methanotrophy by summing all dissolved and particulate methane monooxygenase genes. Finally, we fit the same linear models to individual genes comprising more than 1% of normalized reads of each Subsystem 3 function.

| Patterns in lake chemistry reflect permafrost and landscape position
Lake water chemistry varied both within and among permafrost regions, with within-region variability reflecting local groundwater connectivity. In the sampling period between 2018 and 2019, lakes had relatively high DOC concentrations (17.6-50.5 mg L −1 ), slightly acidic to moderately alkaline pH (6.2-8.5), and moderate to high EC (32-1440 μS cm −1 ) ( Table S1). All the chemical parameters except Fe and Mg concentrations varied among permafrost regions, often by an order of magnitude (Table S2). However, there was at least as much variation among lakes as across permafrost regions for some parameters, including Fe, Ca, S, Mg, Mn, and Zn (Table S2)
We found that the composition of methanogenic archaea varied with the direct and indirect impacts of climate warming, whereas overall microbial community composition, bacteria and non-methanogenic archaea were unchanged. The groupings of climate and hydrological variables were the only ones that explained statistically significant amounts of variation in community composition (Table 1). Together, these two groups explained nearly 50% of the variation of archaea when considering their interactions. The three most common archaeal genera, which collectively accounted for 9%-25% of all archaeal-assigned reads per site, each increased in abundance from the continuous to no permafrost region by a mean ± 95% confidence interval (CI) of 214%-300% (36%-819%; Table S3). All three genera (Methanoregula, Methanosarcina, Methanothrix) also decreased by 99% (90%-100%), 93% (76%-131%) and 96% (70%-100%), respectively, with increasing EC. In support of the interpretation that abundance was limited by the local environment rather than dispersal, we found that the community composition of methanogenic archaea did not vary simply with increasing geographic distance (Mantel's r = .08, p = .168; 0.05% filtering threshold: r = .08, p = .153). The composition of bacteria and all archaea similarly did not vary as sites were more distant from each other (Figure 2c,d).

| Methane cycling genes vary with direct, regional and indirect, local warming effects in permafrost lakes
Consistent with the taxonomic changes, functional genes associated with methane production varied with the direct and indirect impacts of climate warming. The normalized abundance of all genes associated with methanogenesis increased as permafrost cover declined ( Figure 3a). Compared to continuous permafrost, methanogenesis genes were 115% (95% CI: 24%-273%; t 14 = 2.98, p = .010) and 61% (10%-136%; t 14 = 2.67, p = .018) more abundant in the sporadic and no permafrost regions respectively. Genes involved in methanotrophy were similar in all permafrost regions (for all: t 14 < 1.59, p > .135; Figure 3a), consistent with the patterns in bacterial taxonomy ( Figure 2). There was no correlation between the abundance of genes involved in methanogenesis and those involved in methanotrophy, suggesting that methanotrophy did not simply track methanogenesis (ρ = 0.16, p = .51).

F I G U R E 2
Taxonomic composition of microbial communities does not vary with geographic distance. Composition of (a) bacterial and (b) archaeal communities in sediments of 19 lakes across permafrost regions. For archaea, class was reported where available, and phylum (P) otherwise. Actinobacteria and Proteobacteria were dominant bacterial phyla regardless of region. Methanomicrobia and Halobacteria were dominant archaea, with proportional abundance of Methanomicrobia increasing from the continuous to no permafrost region (Table  S2). Pairwise differences in genus-level community composition between lakes for (c) bacteria and (d) archaea calculated with the Horn-Morisita index for genus level. Note the difference in y-axis scale between (c) and (d). We removed all taxa present in a sample with a relative abundance of <0.1%. Results in (c) and (d) were unchanged when communities were filtered at a 0.05% threshold: r = .13, p = .089 and r = .09, p = .134 respectively.
We further disentangled changes in methanogenesis by separately analysing the subset of 64 genes that comprised this functional category in the SEED database. First, we found that the composition of genes varied between permafrost regions and with EC (PERMANOVA: F = 2.99, p = .020 and F = 10.75, p = .002 respectively), which itself explained variation in DOC and ion concentrations that could vary with groundwater connectivity. Sites varied most in composition within the continuous permafrost region suggesting that thaw could be reducing the diversity of functional genes (Figure 3b). Higher EC was also generally associated with communities under sporadic permafrost cover (Figure 3b).

Second, we modelled individual genes that accounted for at least
1% of all the reads associated with methanogenesis. We found that genes encoding both formylmethanofuran dehydrogenase (fmd), unique to methanogenic archaea, and methyl coenzyme M reductase (mcr), which catalyses the final step in methane formation, increased in relative abundance as permafrost cover decreased ( Figure 3c). Other enzymes strongly positively associated with reduced permafrost cover included one-carbon unit carrier methenyltetrahydromethanopterin cyclohydrolase (mch) and the highly salt-dependent formylmethanofuran-tetrahydromethanopterin Nformyltransferase (ftr) (Figure 3c).

| Key functional genes vary with direct, regional and indirect, local warming effects.
In addition to methanogenesis and methanotrophy, we focused on genes associated with organic matter decomposition, sugar utilization and energy production. We identified 27 further Subsystem 3 functions comprising a total of 1798 genes. Although none of these functions varied between the continuous region and either the discontinuous or sporadic region (Table S5) Table S6). Sporadic and no permafrost regions had similar gene composition that was distinct from discontinuous and continuous regions, which also differed from each other (Figure 5b). The regions with lower permafrost cover were more closely associated TA B L E 1 Methanogen genera varied with climate and hydrology. We partitioned the variation in community composition using redundancy analysis. Each partition consisted of multiple indicators of climate: permafrost region, mean annual temperature, mean annual precipitation; hydrology: EC and catchment peatland cover; organic matter (OM) inputs: DOC concentration and SUVA; and lake morphometry (Morph): ratio of lake area to catchment area and maximum lake depth. Values are the variation explained R 2 by each partition considering all two-and three-way interactions with the other partitions, that is, marginal effects, and all partitions together, that is, total model R 2 . Analyses were repeated by discarding genera that comprised <0.10% or <0.05% of normalized counts in a sample. Bolded values are statistically significant at p < .05. For all, n = 19 lake samples. with sucrose phosphorylase (sucP), which is needed to convert sucrose into glycolytic intermediates for microbial utilization, and a TonB-dependent receptor (TDBR) that is predicted to facilitate sucrose uptake (Figure 5b). By contrast, the discontinuous permafrost region clustered with sucrose-6-phosphate hydrolase (scrB) and sucrose operon repressor (scrR) genes, which promote sucrose catabolism, while sucrose utilization genes varied the most within the continuous region (Figure 5b).
We then found that 30 of the 48 most abundant genes across all the Subsystem 3 functions varied between the continuous region and one or more other regions (Figure 5c). Eight of these genes were dehydrogenase enzymes and four were permeases, the latter of which were all more abundant as permafrost cover was lost (Figure 5c).
Likewise, 19 of the 48 most abundant genes were associated with EC in varying directions (Figure 5d). Eleven of these genes were more abundant in lakes with higher presumed groundwater connectivity, including central drivers of sugar metabolism, such as pyruvate kinase (pyk), phosphoenolpyruvate synthase (ppsA) and sugar storage enzymes, such as an archaeal glycogen branching enzyme (glgB) ( Table S8). Finally, although archaea were only a small proportion of the identified reads compared to bacteria, archaea-specific NADPdependent glyceraldehyde-3-phosphate dehydrogenase (gapN) was one of the most abundant genes observed (Figure 5c). gapN is central to the glycolytic pathway, contributing to both glycolysis and gluconeogenesis, and we found that it increased in abundance as permafrost was lost and decreased with increasing EC (Figure 5c).

| DISCUSS ION
Consistent with our predictions, our study suggests climate warming has opposing direct and indirect effects on microbial community structure and function in sediments of permafrost lakes. By sampling 19 lakes across a 1600 km latitudinal transect spanning all major regions of permafrost cover, we assessed the potential impacts of climate warming through a space-for-time substitution. We compared the direct effects of warming that acts regionally, such as on temperatures and the extent of permafrost cover, to its indirect F I G U R E 3 Methanogenesis but not methanotrophy genes are more abundant with decreasing permafrost cover. (a) Partial residuals for the effect of permafrost region on methanogenesis (blue) and methanotrophy (pink) in each of five lakes per region (four in discontinuous permafrost). We separately summed the normalized abundances of all 64 genes comprising the SEED methanogenesis pathway and the two methane monooxygenases that catalyse methane oxidation. Partial residuals visualize differences between permafrost regions while accounting for the effect of groundwater connectivity inferred by EC  (Table S4). Full gene names given in Table S4. Shifts from brown to blue indicate genes that are more abundant as permafrost cover increases. White cells denote functions that do not statistically differ in abundance from the continuous permafrost region (Table S4). Permafrost regions were abbreviated as follows: C, continuous; D, discontinuous; S, sporadic; N, no permafrost.
impacts on local groundwater connectivity and lake water chemistry (Karlsson et al., 2012;Olefeldt & Roulet, 2012;Walvoord & Kurylyk, 2016). We found that metabolically diverse Euryarchaeota methanogens and genes involved in all major methanogenesis pathways responded especially positively to the direct effects of regional climate warming. However, increases in sequences representing methanogens and methanogenesis could be partly offset by greater groundwater connectivity associated with warmer temperatures, as the abundances of these sequences correlated negatively with groundwater connectivity (e.g. Figure 4). These opposing responses to different aspects of climate warming help explain why considerable variation remains unexplained in predictions of methane emissions by Earth system models even when microbial dynamics are considered (Chadburn et al., 2020;Ernakovich et al., 2022). Potential methane production may be further enhanced at warmer sites by the lack of corresponding increase in its oxidation, as methanotrophic gene abundance was relatively consistent across our climate gradient. An increase in methanogenesis is supported by concurrent measurements of methane fluxes at our study sites by Kuhn et al. (2021).
Together, our results provide a new understanding of why methane emissions are highly sensitive to warming temperatures across our wider 260,000 km 2 study region  and elsewhere (Rasilo et al., 2015). Future predictions of methane emissions could now be improved by considering how microbial dynamics vary with the direct and indirect effects of climate warming associated with temperature and hydrology respectively.

| Methanogenesis has contrasting responses to the direct and indirect effects of climate warming
While the direct loss of permafrost cover can create new hydrological connections (Connon et al., 2014), our results suggest that how much these connections impact microbial communities depends indirectly on landscape position ( Figure 6). In regions with lower groundwater connectivity, methanogenesis-inhibiting TEAs such as Fe 3+ , SO 4 2− and NO 3 − are limited (Einsiedl et al., 2008;O'Donnell & Jones Jr, 2006;Pokrovsky et al., 2016). Methanogenesis may therefore benefit from the increasing availability of organic substrates released by permafrost thaw, such as carbon dioxide and acetate ( Figure 6). In support of this interpretation, we found both xylose and sucrose utilization genes, and the most abundant organic matter metabolism genes, were all most abundant in the warmer climate without permafrost. The increasing availability F I G U R E 4 Methane cycling genes are less abundant as groundwater connectivity increases. (a) Electrical conductivity as a proxy for groundwater connectivity was highest lakes in the sporadic (S, n = 5) as opposed to continuous (C, n = 5), discontinuous (D, n = 4) and no permafrost (N, n = 5) regions. Thick grey lines show medians, boxes show interquartile range (IQR), whiskers extend 1.5 times IQR and points are values outside whiskers. (b) Genes involved in methanogenesis (blue) and methanotrophy (pink) decreased with electrical conductivity across 19 lakes. We summed all 64 genes comprising the SEED methanogenesis pathway and the two methane monooxygenases that catalyse methane oxidation. Lines are mean effect of electrical conductivity ±95% confidence intervals. (c) The most abundant methanogenesis genes generally decreased with conductivity across the 19 lakes. White cells indicate effects that are not statistically significant (Table S4). Gene abbreviations and model fits as in Figure 3.
and digestion of organic matter previously stored in permafrost could further explain our results in the sporadic and discontinuous permafrost region (Feng et al., 2020). In the non-permafrost region, higher within-lake productivity (Karlsson et al., 2009) and fresh peat-derived inputs could additionally stimulate decomposition (Burd et al., 2020).
By contrast, methanogenesis may be more strongly inhibited by TEAs in regions with higher groundwater connectivity ( Figure 6).
As groundwater connectivity (measured by EC) can also promote sugar metabolism genes that produce substrates for methanogenesis (Costello et al., 1991), it may have also mixed and weaker effects on methanogenesis than temperature. EC was also only higher in the sporadic permafrost region. This result possibly reflected the location of our study lakes in a groundwater discharge zone given their proximity to major rivers and surrounding areas of higher elevation. Together, the magnitude of the predicted net increase in methanogenesis under future climate warming , will be influenced by local surficial geology and associated groundwater characteristics ( Figure 6).

| Methanotrophy only indirectly varies with climate warming
Our results suggest methanogenesis will become increasingly dominant over methanotrophy as northern lakes warm ( Figure 6).
The abundance of microbial genes involved in methanotrophy was negatively associated with groundwater connectivity but not regional climate warming. As methanotrophy is highly dependent on methane concentrations (Mayr et al., 2020), it may have simply decreased as methanogenesis was inhibited by the increased delivery of TEAs from mineral groundwater. Some methane oxidation will also be obligate anaerobic and proceed through reverse methanogenesis in certain Euryarchaeota via the enzymes encoded by mcr F I G U R E 5 Genes involved in organic matter metabolism vary between permafrost regions and with groundwater connectivity. (a) Subsystem 3 level functions involved in organic matter metabolism, for which normalized counts varied between lakes in continuous (C, n = 5) and no permafrost (N, n = 5) regions. These functions did not vary between the continuous and any other region (Table S5). Thick grey lines show medians, boxes show interquartile range (IQR), whiskers extend 1.5 times IQR and points are values outside whiskers. (b) Nonmetric multidimensional scaling (NMDS) of 33 genes involved in sucrose utilization in each of 19 lakes denoted by points (Table S7). Coloured lines connect lakes in the same region. For visualization purposes, we only plotted the 10 most abundant genes within the sucrose utilization subsystem in red. NMDS stress = 0.06. For organic matter metabolism genes that accounted for ≥1% of reads within their respective Subsystems, we plotted those that (c) differed in one or more regions relative to continuous permafrost and (d) varied positively (green) and negatively (brown) with electrical conductivity as a proxy of groundwater connectivity across 19 lakes. Effects were estimated with linear models with corresponding R 2 values of the overall model in parentheses. Only statistically significant effects are coloured. Unabbreviated gene names and model outputs are in Table S8.
and mer (Evans et al., 2019), rather than methane monooxygenases in bacteria (Haroon et al., 2013;Timmers et al., 2017). However, had substrate availability and reverse methanogenesis alone been important, methanotrophy should have increased from continuous to no permafrost cover, as methane substrates and mcr and mer genes became more abundant. The absence of this response in abundances of methanotroph genes suggests that methane oxidation was primarily performed by bacteria under aerobic conditions rather than anaerobically by archaea (Deutzmann & Schink, 2011;Eller et al., 2005).

| Organic matter processing genes have mixed responses to climate warming
Metabolic flexibility likely contributed to high levels of functional redundancy and little change in bacterial taxonomic composition with climate warming or geography (Ernakovich et al., 2022). The most abundant organic matter metabolism genes, such as xynA (endo 1,3-beta-xylanase A precursor), paaC (3-hydroxyacyl CoA dehydrogenase) and rbsA (ribose ABC transporter ATP-binding protein), many of which encode monosaccharide-degrading enzymes, generally increased as permafrost cover was lost and more labile organic matter became available (Dutta et al., 2006).
However, there was no consistent change in these genes with groundwater connectivity, and most genes associated with organic matter metabolism and energy production also did not vary with permafrost cover. Bacteria that harbour most of these nonmethanogenesis genes likely have mixed responses to the different electron acceptors delivered from groundwater. For example, SO 4 2− reduction can be inhibited by non-O 2 TEAs such as Fe 3+ (Lovley, 1991), and other bacteria can use TEAs for carbon metabolism (Chen, Leung, et al., 2021). Several genes that we found were positively associated with connectivity also encode enzymes involved in sugar catabolism that rely on metal cofactors delivered by groundwater, for example, adhE (bacterial alcohol dehydrogenase), ppsA (phosphoenolpyruvate synthase) and pykF (pyruvate kinase) (Steinmann & Shotyk, 1997). Increasing groundwater connectivity may also deliver and enhance within-lake production of bioavailable organic matter (Olefeldt et al., 2013), explaining the positive response of many sugar utilization genes. Therefore, changes in taxonomic composition may only arise if the functions selected by the environment are strongly phylogenetically conserved, as we found with methanogenesis and archaea. More Bioinformatics processing of sequences was performed by Lucas P.
P. Braga. Statistical analysis was performed by Johanna C. Winder

F I G U R E 6
Climate warming has opposing direct and indirect effects on microbial processes associated with methane emissions. Left: In the absence of warming, intact permafrost stores organic matter (OM) and acts as a barrier to deep groundwater flows that deliver terminal electron acceptors (TEAs), especially to lower elevations. Right: Permafrost thaw arising from climate warming makes OM available for decomposition, which we found can shift microbial communities towards methane production. Loss of permafrost may offset some of the expected shift towards methanogenesis at lower elevations by removing barriers to the delivery of TEAs from upstream mineral soils, but also reduce methane oxidation. and Andrew J. Tanentzap. Johanna C. Winder wrote the paper with contributions from Andrew J. Tanentzap. All authors commented on the manuscript.

ACK N O WLE D G E M ENTS
We

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no competing interests relevant to this study.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in multiple locations. Metagenomic data generated during this study are openly available in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB36558 (https://www.ebi.ac.uk/ ena/brows er/view/PRJEB 36558). Seasonal and daily data including lake coordinates and chemical measurements are openly available in Table S1 and in dryad at https://doi.org/10.5061/dryad.sj3tx 968t. Estimates of peatland lake area were extracted from the Boreal Arctic Wetland and Lake Database (BAWLD) and are openly available at the Arctic Data Center (https://doi.org/10.18739/ A2C82 4F9X).