Spatial heterogeneity of methane concentration, diffusion and production potential, and methanogen community abundance and diversity in sediments from Lake Remoray (France)

Lake ecosystems contribute signicantly to atmospheric methane and are likely to become even bigger methane emitters with the global spread of hypoxia/anoxia in freshwater ecosystems. Here we characterized the spatial heterogeneity of methane production potential, methane concentration, archaeal and bacterial and methanogen communities across Lake Remoray sediment during the summer period, when hypoxic conditions settle in the deepest part of the water column. It was hypothesized that methane concentration and production would be higher in the deeper part of lake, our results showed that some littoral areas exhibited similar or higher values than the deepest area. The full 16S rRNA gene sequencing dataset counted 41 methanogen OTUs in abundances that depended more on sampling-site location than on the water depth gradient. The methanogen co-occurrence network revealed the existence of 5 distinct sub communities, suggesting the presence of different methanogen niches across Lake Remoray. The variation in abundance of the two larger methanogen sub communities was signicantly related to methanogenesis potential and sediment methane concentration across-lake. Taken together, our results highlight the need to consider the ecosystem at whole-lake scale when calculating methane budgets, and the need to couple spatial and temporal approaches in order to better estimate lake ecosystems contribution to atmospheric methane.


Introduction
Methane ranks just behind carbon dioxide as a major greenhouse gas, with an 12-year lifetime in the atmosphere, a global warming potential 72 times higher than CO 2 over a 20-year timescale (25 times higher over 100 years), and a radiative forcing of + 0.48 W m -2 [1]. Methane concentration in the atmosphere is currently about 1850 ppb [2], which is an increase of more than 150% since the preindustrial era and rising fast at the unprecedented rate of > 5 ppb yr -1 in the 2004-2017 period [2].
Atmospheric methane is mostly anthropogenic: fossil fuel extraction and agricultural practices account for about 73% of methane emissions. Other sources are natural: wetlands ecosystems are the biggest natural source of atmospheric methane, followed by lake ecosystems which account for 6-16% of nonanthropogenic emissions [3], i.e. up to 76 Tg C-CH 4 emitted to the atmosphere per year [4]. Recent models project that lake emissions will increase under the combined effects of rising temperature [5,6] and global-scale eutrophication of water bodies [7]. The isotopic signature of atmospheric methane is also changing, as δ 13 C-CH 4 has fallen from -47.06‰ to -47.36‰ in less than 10 years, likely re ecting changes in both sources and sinks of atmospheric methane [2].
Microbial actors of the methane cycle are the methanogenic archaea, which are involved under anoxic conditions in the nal stage of reducing organic matter (OM) into methane [8], and methanotrophic organisms (archaea and bacteria) that oxidize, with e ciencies ranging from 30% to 99% in lake ecosystems [9],) methane under both anaerobic and aerobic conditions into carbon dioxide [10,11]. In lakes, methanogenesis mainly takes place in the sediment compartment and in the deeper water column layers under temporary or permanent anoxic conditions [12]. Methane oxidation occurs aerobically when carried out by methanotrophic bacteria or anaerobically in the presence of sulfates or nitrates as terminal electron acceptors when carried out by methanotrophic Archaea [13]. If not oxidized, methane can be emitted to the atmosphere through four well-described pathways: bubbling, diffusion, by lake overturn when methane trapped in the hypolimnion is released, and via ows related to aquatic plants and roots that facilitate sediment-water column gas exchanges [3].
Among the abiotic factors likely to affect the methane cycle, those that vary spatially and seasonally and with strati cation of the water column will de ne the location of the oxic/anoxic interface and in uence methanogenic and methanotrophic community diversity and activity rates [14][15][16][17]. Methane oxidation dynamics are in uenced by concentration of methane [18], oxygen [19][20][21] or nitrogen [21] and by methanotroph grazing by benthic invertebrates [22], whereas methanogenesis depends more on redox potential [23], quality and quantity of available OM [9], and temperature [20]. In particular, it has been shown that methane production in lake sediments is correlated with the C/N ratio of OM: methane production rates are higher at C/N < 10 [20]. This suggests that methanogenic and methanotrophic communities are likely to be sensitive to variations in redox potential but also to qualitative and quantitative variations in OM inputs. At lake scale, seasonally-driven water column strati cation maximizes the spatial occurrence of favourable redox methanogenesis conditions in the hypolimnion [24], and thermal strati cation is expected to last for longer periods in future global climate change scenarios [25][26][27].
Here we aimed to assess the spatial heterogeneity of methane concentrations, diffusive uxes and potential production rates in the sediment compartment of Lake Remoray (France). Lake Remoray is characterized by eutrophication due to OM, nutrient and contaminant inputs from the watershed, and by long periods of hypoxia in the deep water column (Belle et al., 2016). In July 2017, sediment cores were collected from 10 sampling sites across the 95-ha lake. Environmental descriptors were related to physicochemistry of the sediment compartment (OM content). Abundances of archaeal and bacterial communities involved in methane production and oxidation were assessed by targeting mcrA (methyl coenzyme M reductase) and pmoA (particulate methane monooxygenase) as functional markers. Diversity, structure and composition of sediment archaeal and bacterial communities were assessed using an Illumina sequencing approach targeting 16S ribosomal RNA (rRNA) genes. The whole sequencing dataset served to characterize the methanogen communities by co-occurrence network analysis.
We expected that (i) methane concentration, production and diffusion would be higher in sediment from the deeper zones due to extensive recent OM accumulation and to the longer time period of anoxia in overlying water column layers, and that (ii) spatial variations in sediment methane concentration, production and diffusion would be related to the diversity and to their metabolic genetic potentials (i.e. abundance) of the in-sediment methanogenic communities.

Materials And Methods
Study site and core collection Lake Remoray (46°46′12″ N; 6°15′49″ E; Fig. 1) is located in the Jura Mountains (eastern France). The water depth of this moderately-sized lake (95 ha) reaches 27 m. The lake is currently meso-eutrophic, with an accumulation of organic matter in the central and deepest part of the lake [29]. It is fed by two tributaries (the Drésine and the Lhaut), has one outlet (the Taverne) (Fig. 1), and is bordered by a farm and a campsite. The catchment basin of Lake Remoray consists of 46.7% forest (mainly coniferous forest), 42.5% agricultural land, 5.5% peat bog, 2.6% water area, and 2.7% urban area [28].
In 2017 (July 4 th and 5 th ), short 90-mm-diameter sediment cores were collected at 10 sampling sites located all over the lake using a UWITEC gravity corer (Fig. 1) (Fig. 2). At sampling time, water temperature at the water -sediment interface was of 5°C for the deeper sites and 21.7 °C for the shallower sites. Two littoral sites were included in the tributaries area to include the potential in uence of the Drésine and the Lhaut on the lake. Two cores were collected at each site, and the rst 10 centimetres of sediment from the sediment water interface in each core was sub-sampled directly in the eld. For the rst core, sediment samples were collected through a pre-pierced liner using precut autoclaved syringes at determined sampling depths (n=10, every centimetre for the rst 10 centimetres), and used for biogeochemical measurements (CH 4 concentrations and diffusive ux calculations) and microbiological analyses (molecular-based analyses). For the second sediment core, sediment samples were extruded using a core extruder and 2-cm sediment sections were collected using a core cutter (n=5, every two centimetres for the 10 rst centimetres) to serve for biogeochemical measurements (potential activity rates) and sedimentological analysis.

Sediment characteristics
Sediment particle size distribution was determined by laser diffraction on dry sediments (Mastersizer S, Malvern, United Kingdom), and particles were grouped by median grain size. Sediment water content was estimated by drying at 60°C for 48 h. Organic matter (OM) content was measured by loss-on-ignition [30].
Organic carbon (OC) content and porosity were calculated from water and OM contents using formulas provided elsewhere [31].

Methane concentrations and uxes
Sediment samples (0.5 mL) subsampled from the sediment core every centimetre through pre-pierced holes were placed into glass vials with NaOH (4 mL of 2.5% NaOH in 20-mL vials), covered with a butyl stopper, and sealed with an aluminium crimp. Dissolved CH 4 concentrations were measured in the headspace by gas chromatography using an Agilent 490 MicroGC thermal conductivity detector equipped with a PoraPLOT Q column with helium as vector gas.
Methane uxes (Js) were calculated based on Fick's rst law using the top 10 cm of measured CH 4 concentrations. The equation from Maerki et al., (2004) was used to correct for porosity and tortuosity. The CH 4 diffusion coe cient, D, was interpolated by sampling site water temperature at the watersediment interface according to Lerman (1979). Js was calculated according to Sollberger et al., (2014) as follows: Methane production potential CH 4 production potentials were measured as described in Fuchs et al. (2016) by incubating lake sediment under controlled conditions. Immediately following eld sampling, wet sediment from each targeted 2-cm section of extruded sediment was placed in a ask, lled to the top to exclude air capture, sealed, and stored at 4°C until potential activity measurements. For CH 4 production potential, 5-10 g of wet sediment was transferred to 150 mL asks added with 10 mL of demineralized water. The asks were hermetically sealed with rubber caps, inerted with helium, and incubated at 20°C under constant moderate agitation during the course of the experiment. Methane concentrations were then measured by gas chromatography (n=1) at t0 = 1 h (allowing 1 h of medium equilibration) and at 5 h using the Agilent 490 MicroGC thermal conductivity detector. For each ask, dry weights of incubated sediments were determined after drying at 105°C for 48 h. Methane production was expressed as ng of C-CH 4 g −1 sediment dry weight [DW] h −1 based on the increase in CH 4 concentrations.

Molecular analyses
Samples for nucleic acid-based analyses were collected using pre-cut autoclaved syringes, transferred to cryovials, immediately stored on liquid nitrogen in the eld, taken back to the lab and held at -20°C until nucleic acid extraction was realized. Total DNA was extracted from 0.5 g of wet sediment using a NucleoSpin Soil Kit following the manufacturer's instructions, and using SL1 lysis buffer and additive Enhancer SX buffer (Macherey-Nagel, Hoerdt, France). The extracted DNA was quanti ed uorometrically after staining with QuantiFluor dsDNA Dye (QuantiFluor dsDNA System, Promega, Charbonnières-les-Bains, France) using a Plate Chameleon™ uorometer (Hidex, Turku, Finland; excitation: 485 nm, emission: 590 nm).
The abundances of total and functional organisms were assessed by quantitative PCR targeting the 16S rRNA, mcrA and pmoA genes. initial cycle of 95°C for 10 min, 45 cycles of 95°C for 15 sec, 60°C for 30 sec and 72°C for 30 sec, and then a melt-curve step (55°C to 95°C). A plasmid containing a single copy of 16S rRNA genes ampli ed from Escherichia coli (Bacteria) was diluted from 10 8 to 10 1 copies per assay and used in triplicate to produce the standard curve. Ampli cation e ciency was 105%. Methanogen and methanotroph abundance was evaluated by quantitative PCR on mcrA and pmoA genes using the protocol described in Fuchs et al. (2016) but with a nal qPCR mix adjusted to 25 μL. Reactions were performed in a nal volume of 25 μL containing 1× Brilliant II SYBR Green QPCR Master Mix (Agilent), 0.3 mg mL −1 bovine serum albumin (Sigma-Aldrich), 0.5 μM of primers ME1 and ME2 [36] and A189f and mb661r [37] for mcrA and pmoA assays, respectively, and 2 ng of template DNA. A plasmid containing a single copy of mcrA gene ampli ed from a Methanosarcina thermophila a liated clone (KR011363) and pmoA gene ampli ed from Methylomonas methanica DSMZ 25384 DNA extract was diluted from 10 7 to 10 1 molecules per assay. Samples and standards were used in triplicate. Ampli cation e ciencies were 93% for mcrA and 99% for pmoA. One or two copies of mcrA per genome is reported [38] while two copies of the pmoA gene were assumed to be present per cell [39].
PCR ampli cation for high-throughput 16S rRNA sequencing was carried out with the universal primer pair 515F (5′-GTGYCAGCMGCCGCGGTA-3′) and 909R (5′-CCCCGYCAATTCMTTTRAGT-3′) targeting the hypervariable V4-V5 region of the 16S rRNA gene [40]. Indexes were integrated to both primers following the dual-indexing procedure described by Kozich et al. (2013). Triplicate PCR ampli cations for each sample were realized with a total amount of ~5 ng of DNA per reaction. Amplicons were quanti ed using a PicoGreen™ assay (Life Technologies, Carlsbad, CA) and pooled equimolarly. The nal pool was puri ed using CleanPCR beads (CleanNA). Sequencing was realized by Fasteris (Geneva, Switzerland) on an Illumina HiSeq 2500 system with 2 × 250 bp. Analysis yielded 6.7 Gb of sequences with average error rate of 0.822% and an average Q30 of 90.3%.

Bioinformatics analysis
Adapter sequences were removed using Trimmomatic [42] and barcode-sorted using a Fasteris internal script (Perl script). Sequences were then processed using the FROGS (Find Rapidly OTUs with Galaxy Solution) Galaxy-supported pipeline [43]. Paired-end reads were joined using FLASH [44] and qualitychecked using FastQC. Sequences with non-mismatching primers were kept and then ltered by size (between 350 and 500 bp), and those containing N bases were discarded. The 16S rRNA gene sequences were then denoised and clustered using the Swarm method [45] with a maximum 3-bases difference. Clusters having less than 0.005% abundance and with occurrence in less than two samples of the total library were deleted. Chimeras were removed using the vchime tool of the vsearch package [46]. A liation was done using the Silva SSU database v.123 [47] via BLAST [48] with multiple a liation allowed and manual curation. All analyses were done on the Galaxy instance of the INRA MIGALE bioinformatics platform (http://migale.jouy.inra.fr). The raw datasets are available in the European Bioinformatics Institute (EBI) database system (in the European Nucleotide Archive) under project accession number PRJEB 43875.

Statistical analysis
Methane production potential, methane concentration, methane diffusion, microbial density, mcrA/16S rRNA genes ratio and pmoA/16S rRNA genes ratio were compared between the littoral and profundal zones of Lake Remoray using a nonparametric Kruskal-Wallis test for each stratum of sediment core samples.
Based on the taxonomic classi cation, methanogen sequences were extracted from sequencing dataset.
We used all the orders of methanogen known in literature and which represented by Methanobacteriales, Methanomassiliicoccales, Methanomicrobiales, Methanosarcinales, Methanomethyliales and Methanofastidiosales [8]. Spatial variations in archaeal and bacterial, and methanogen community structure across Lake Remoray were characterized using UniFrac distance [49]. Nonmetric multidimensional scaling (NMDS) was used to graphically depict spatial variations between the archaeal and bacterial, and methanogen communities. The signi cance of the observed clustering of samples on the ordination plot was assessed by analysis of variance using distance matrices (ADONIS test, 999 permutations).
We used the 'DESeq2' R package [50] to calculate the signi cant changes in methanogen composition (log2 fold change in relative abundance of each OTU) between the littoral and profundal zones of Lake Remoray. Differences in methanogen OTU abundances between in-lake depth zones were considered signi cant if their adjusted p-values were <0.01.
The methanogen community co-occurrence network was built using the 'igraph' R package. Pearson correlations between all pairs of 41 OTUs (OTUs found in at least 15 samples) were calculated, and the pvalues were adjusted using Benjamini-Hochberg (BH) method [51]. Only signi cant correlations (p<0.05, post-BH correction) with R>0.5 were graphed. Network modules (subcommunities in the network) were identi ed using the 'WGCNA' R package [52] on the adjacency matrix of Pearson correlation values between each pair of OTUs with a thresholding power of 6. The adjacency matrix was transformed into a topological overlap matrix (TOM) to minimize possible errors, and the OTUs were clustered in modules with a minimum number of three OTUs per module. Modules having <0.3 dissimilarity in eigengenes were subsequently merged, and Pearson correlations between merged modules and the methane functional parameter dataset (CH 4 production, CH 4 and OM concentration, and mcrA/16S rRNA genes ratio) were calculated.
All statistical analyses were performed using with the free R software (version 3.5.3).

Sediments characteristics
Sediment texture was relatively similar between the different sampling sites of Lake Remoray, with about 73% silt, 25% sand and 2% clay (Fig. 2). Only sampling site #1413 close to the Lhaut tributary had a different sediment texture, with 19% silt, 81% sand and less than 1% clay.
Organic matter (OM) and organic carbon (OC) content in the rst 10 centimetres of sediment cores varied from 6.3% to 38% and from 3.7% to 22%, respectively, depending on sampling-site location (Fig. 2). The highest OM and OC contents were found at sampling site #1414 close to lake tributaries, with in average for the rst 10 cm of the sediment core containing an average 24% OM and 14% OC. The deepest sampling site (#1441) was characterized by 12% OM content and 7.8% OC content.
Methane production potential, concentration and diffusion across Lake Remoray Methane production potential, concentration and diffusion were not related to water depth at the sampling site (littoral vs profundal zone) and were not signi cantly different between stratum depths along sediment cores (Kruskal-Wallis test in all sediment strata; data not shown) (Fig. 3).
In the profundal zone, methane production potential was maximal in the rst 2 cm of sediment cores (Fig.   3a). The highest methane potential production was recorded at the deepest sampling site, with 6563 ng C-CH 4 g DW -1 h -1 . Methane production potential decreased gradually along the strata of sediment cores at all profundal-zone sampling sites. Average methane production potential in the littoral zone varied from 789 to 1674 C-CH 4 g DW -1 h -1 along the sediment core strata. The highest littoral-zone methane production potential was recorded at sampling site #1419, with 1833 to 3263 C-CH 4 g DW -1 h -1 over the rst 10 cm of the sediment core. The lowest methane production potential was recorded at sampling sites #1413 and #1414, with an average of 179 and 705 C-CH 4 g DW -1 h -1 , respectively, over the rst 10 cm of the sediment core.
Methane concentrations increased with sediment core depth (Fig. 3b), ranging from 2 to 8.2 mg CH 4 L -1 for littoral-zone cores and 2.6 to 11 mg CH 4 L -1 for profundal-zone cores. The deepest sampling site (#1441) had the highest methane concentrations along the sediment core. Sampling site #1413, in the tributaries area, had the lowest methane concentrations in the sediment core.
Genetic potential of methane cycle functioning across Lake Remoray Whatever the sampling-depth zone, bacterial densities assessed as 16S rRNA gene copies g -1 DW (Fig.  4a) decreased from top to bottom of the sediment cores. No signi cant differences were observed between littoral and profundal zones (con rmed by Kruskal-Wallis tests for all sediment strata) and the variation of bacterial densities depended more on sampling-site location across Lake Remoray.
Gene-based potential methanogenesis (mcrA/16S rRNA genes ratio) and methanotrophy (pmoA/16S rRNA genes ratio) (Fig. 4b and 4c) varied depending on sampling-site location across Lake Remoray. There were no signi cant differences between littoral and profundal zones, except pmoA/16S rRNA genes ratio which was signi cantly higher in the 2-cm stratum of littoral-zone sediment cores. The microbial community had 700-fold higher methanogenesis potential than methanotrophy potential. The highest methanogenesis potential was found for profundal-zone site #1440 (average mcrA/16 ratio of 0.31 over the whole sediment core) and for littoral-zone sampling site #1398 (average mcrA/16 ratio of 0.34 over the whole sediment core). The lowest methanogenesis potential was found for littoral-zone sampling sites #1413 and #1414 (average mcrA/16S rRNA genes ratio of 0.09 and 0.08, respectively, over the whole sediment core).
Archaeal and bacterial, and methanogen community structure across Lake Remoray The high-throughput 16S rRNA gene sequencing approach applied here generated 3,788,355 quality sequences and recovered 1566 OTUs a liated to both the Archaea and Bacteria domains.
NMDS analysis of the full archaeal and bacterial sequence dataset highlighted a distinct difference in sediment community structure between littoral and profundal zones (Fig. 5a) (ADONIS test, R²=0.13, P=0.001). This difference was associated with different taxonomic compositions: Proteobacteria, Bacteroidetes, Planctomycetes, Acidobacteria and Chloro exi were more abundant in the profundal zone whereas Cyanobacteria, Spirochaetes and Euryarchaeota were more abundant in the littoral zone of Lake Remoray. There was no effect of sampling depth along the sediment cores on microbial structure (ADONIS tests; R²=0.10, P=0.911).
Based on the taxonomic classi cation in the 16S rRNA gene sequences dataset, we identi ed 41 OTUs a liated to the methanogens group, representing 6.11% of the total sequences. NMDS analysis of methanogen-sequences dataset revealed two distinct groups on the ordination plot (Fig. 5b). The rst group was represented by sediment cores collected in the tributaries area, and the second group clustered the sediment cores from all other sampling sites across the lake. ADONIS tests showed a signi cant effect of sampling depth between littoral and profundal zones (R²=0.22, P=0.001) and no signi cant effect of sampling depth along the sediment cores (R²=0.03, P=0.999) on methanogen community structure.
DESeq2 analysis on the 41 methanogen OTUs identi ed methanogen genera accounting for the differences between littoral and profundal zones (Fig. 6). Profundal-zone sediments had a higher relative abundance of OTUs a liated to Methanomicrobiales order (81% of methanogen sequences) compared to littoral zone (65% of methanogen sequences). Sediments in the littoral zone had a two-time higher relative abundance of OTUs a liated to Candidatus Methanoperedens, C. Methanomethylicus and the Methanosaeta genus.

Linking methanogen communities and CH 4 -cycling parameters across Lake Remoray
The methanogen community co-occurrence network was built based on pairwise correlations between the 41 OTUs a liated to methanogens. The co-occurrence pattern within the methanogen community was determined by only keeping signi cant Pearson's correlations. The resulting network ( Fig. 7a and b) was composed of 131 edges (all edges were positive correlations) representing 15% of possible withinnetwork connections (network density). The topological network parameters were: average degree=6.4, diameter=6, modularity=0.48, transitivity=0.73, average path length = 2.46.
Covariation similarity analysis between methanogen OTUs showed 5 distinct modules within methanogen network and containing more than three OTUs (Fig. 7b). Analysis of the proportion of methanogen OTUs composing the network modules of Lake Remoray (Fig. 7c) showed that the methanogen OTUs of modules 1 and 2 were dominant in sediments close to the lake tributaries. The methanogens of module 5 dominated the profundal zones while the OTUs of module 3 were more abundant in littoral zones. The methanogens composing module 4 were found speci cally in profundal zones, especially in the deepest sediment of Lake Remoray.
Analysis of module-functional parameter relationships based on the 100 samples collected from the 10 sample cores (Fig. 7d) revealed that sub-communities of modules 1 and 2 were only signi cantly correlated to OM concentration, whereas modules 3 and 4 (representing the largest methanogen subcommunities) were negatively correlated to OM and signi cantly positively correlated to CH 4 concentration and functional methanogenesis potential. Methane production in Lake Remoray was signi cantly correlated to the abundance of module-3 and module-4 sub-communities.

Discussion
Lake Remoray presents long periods of hypoxia in the deepest part of the water column [28] due to an excess of autochthonous material deposition driven by primary overproduction and, to a lesser extent, to organic material inputs carried by its two tributaries (Tardy et al., In press). Methanogen species a liated to the Methanomicrobia class were shown to be widely distributed in the sur cial sediment layers of lake Remoray deepest parts (Tardy et al., In press), allowing to expect methane production in the sediment and methane diffusion along the water column, at least during periods of hypoxia. The aim of this study was to characterize the spatial variation of methane production across Lake Remoray sediment during the summer period, when dissolved oxygen levels decrease sharply with depth in the water column. Due to seasonally anoxic periods in profundal sediment and differences in OM content and quality, our main hypotheses were that it would select for speci c and more abundant methanogen communities within the sediment microbial assemblages, and that methane production potential would increase from the littoral zone to the deepest zone of Lake Remoray.
Spatial heterogeneity of methane production and concentration across Lake Remoray Both methane production potential and methane concentration were more related to sampling-site location than to water depth of the collected sediment cores. Similar or even higher potential methane production rates were recorded in sediment collected in the littoral zone (close to the farm and outlet) than in sediments collected in the deeper zones, with the exception of littoral-zone sampling sites located in the tributary plumes that had the lowest values for methane production potential.
The literature shows that methane production rates depend on the amount, nature and composition of OM inputs [20,[53][54][55]. In a controlled experiment, West et al. (2012) found higher methane production rates following algal biomass input than terrestrial carbon addition whereas methanogenesis potential was similar between conditions. Duc et al. (2010) studied sediments of 8 boreal and temperate northern lakes and found higher methane production rates in sediments characterized by a C:N ratio below 10, whatever the different incubation temperatures tested. In Lake Remoray, sediment had a higher C:N ratio (average C:N ratio=12) in the tributary area than in any other location of the lake (average C:N ratio=8) (Tardy et al., In press).
As previously observed [57,58], methane production potential was higher in the top layers of the sediment cores and ranged from 146 to 6563 ng C-CH 4 g DW -1 h -1 in the rst 6 cm of the core (Fig. 3).
Methane production potential measured in Lake Remoray spanned a higher range than rates measured in Lake Onego (0-456 ng C-CH 4 Duc et al., 2010). Even though the incubation times used here were short (5 hrs), methane production is still temperature-dependent [5,59,60] and this could have in uenced the values obtained. However, the incubation temperature used here was chosen to re ect the temperature recorded at the sampling site in the shallower areas (21.7°C) of the lake.
The highest potential methane production rates were in the sediment top layers. However, the methane concentrations measured in situ in the sediment core layers increased with distance from the watersediment interface at every sampling location. This could be related to (i) sub-favourable in situ conditions for methanogenesis, likely related to oxygen concentration or redox potential, and/or to (ii) methane loss from sediment to water column though either oxidation or emission, or both. Conceptual models based on methane concentration vs. depth curves could help to estimate the balance between oxidation and emission [61], but they are only qualitative and remain di cult to interpret. The analysis of the carbon isotopic signature of methane would have helped to quantify the contribution of CH 4 oxidation [19] to diffusion mitigation. Methane oxidation can occur under both aerobic and anaerobic conditions, and is a metabolic pathway supported by both archaeal and bacterial microbial species [12]. Here, aerobic methanotrophy community abundance, assessed by quantifying relative abundance of the pmoA gene, represented only 0.01% to 0.03 % of the sediment bacterial community. Although only sparsely represented, aerobic methanotrophs have very e cient enzymatic activities and can oxidize signi cant amounts of methane even at low populations densities (Mayr et al. 2020a;2020b;Reis et al., 2020). However, as oxygen diffusion in the sediment compartment is expected to be low, aerobic methane oxidation may not be the only pathway involved in Lake Remoray, especially in the deeper area of the lake and/or deeper layers of the sediment column, even though methanotrophs are also active in anoxic conditions [65]. Anaerobic oxidation of methane is coupled with nitrate, sulphate or iron reduction [12]. The occurrence of the anaerobic methanotrophs Ca. Methanoperedens-related sequences suggests that the anaerobic oxidation of methane using nitrate as terminal electron acceptor via the reverse methanogenesis pathway [66] could take place in Lake Remoray, as previously described in other lakes [17,67]. Further studies, based both on pore-water chemistry and carbon isotopic signature of CH 4 , are needed to better estimate how different oxidation pathways contribute to methane emission mitigation.
In Lake Remoray, methane diffusion from the sediment compartment to the water column varied Taken together, our results argue against our main hypothesis and also highlight the signi cant contribution of littoral areas to methane production in this lake system. Previous studies have reported higher methane production in littoral zones relative to profundal-zone sediment [70][71][72]. This nding has important implications, since global estimations of methane emission from lake ecosystems are mainly based on a single measurement carried out in the deep areas of lakes [9,[73][74][75].
Spatial variation of sediment archaeal and bacterial, and methanogen communities As observed in other freshwater lake ecosystems, water depth signi cantly affected sediment archaeal and bacterial community composition [76][77][78][79]. Littoral-zone sediment was dominated by Cyanobacteria species (29.7% of total sequences) while profundal-zone sediment was dominated by Proteobacteria, Bacteroidetes, and Chloro exi (25.4, 21.2 and 10.5% of total sequences, respectively). There was no signi cant effect of in-core stratum depth on archaeal and bacterial structure and composition. This nding contrasts with numerous studies showing a vertical structuration of bacterial and archaeal communities with the sediment depth in lake ecosystems [80][81][82][83][84][85]. In these studies, archaeal and bacterial composition changes were associated with different patterns of variation in chemical properties along the sediment depth, such as OM quality, oxygen concentration or redox potential. The lack of vertical structuration of sediment archaeal and bacterial communities found here suggests that chemical properties in the rst 10 centimetres of the sediment compartment are relatively homogeneous.
Based on the taxonomic a liation in the 16S rRNA gene sequences dataset, we identi ed 41 OTUs a liated to archaea involved in the methanogenesis function. The relative abundance of methanogen sequences varied from 1 to 17% according to location and in-core sediment depth. A large proportion of sequences was a liated to Methanoregula (67% of total methanogen sequences). Methanoregula is a hydrogenotrophic genus able to reduce CO 2 with hydrogen as electron donor, and it has already been reported as predominant in many freshwater lakes (Borrel et , 2020). The second most abundant methanogen was Methanosaeta (15% of total methanogen sequences), an acetotrophic genus able to produce methane with acetate as nal electron donor [88]. Based on the DNA abundance ratio of these two methanogenic populations in sediment, it suggests that methanogenesis in Lake Remoray results mainly from the hydrogenotrophic pathway and to a lesser extent the acetoclastic pathway, but this would require to be con rmed by activity measurements. A previous study suggested that the hydrogenotrophic pathway outcompetes the acetoclastic pathway in lake sediment if OM is recalcitrant [58].
Methanosaeta was the more abundant methanogen in the littoral zone whereas species a liated to Methanomicrobiales and Methanofastidiosales were the more abundant methanogen in profundal sediments. Littoral sediments located in the tributary area had a speci c methanogen composition. In this part of lake, the methanogen community was largely associated to OTUs a liated to Methanobacterium and Methanosarcina, highlighting the signi cant horizontal structuration of the methanogen community at the whole-lake scale. A few studies have characterized the horizontal distribution of sediment methanogen communities within a single lake ecosystem. Only Berberich et al.
(2020) in Lake Harsha showed large variation in methanogen community composition and structure between riverine, lacustrine and transitional areas of the lake, as well as a higher proportion of Methanosarcina species in the riverine zone.
To go further in the characterization of the methanogen communities, a co-occurrence network was built from the 41 OTUs detected in all sediment samples. The resulting methanogenic network highlighted exclusively positive relationships between the different OTUs, suggesting the absence of competition and the predominance of facilitative interactions within the sediment methanogen community of Lake Remoray. The other interesting metric from the methanogen network is modularity, which was positive and high enough (Q=0.48) to point to the possible presence of community structures (or 'modules') within the network [89]. In biotic co-occurrence networks, the presence of sub communities is usually related to the heterogeneity of habitat and to the diversi cation of ecological niches in the studied ecosystem [90,91]. Our analysis in Lake Remoray sediment found 5 distinct methanogen sub communities (modules) in which interactions are more frequent than with the rest of the community. Two of these sub communities were particularly large (modules 3 and 5; Fig. 7) and both grouped the interactions between OTUs a liated to Methanoregula, Methanolinea, Methanosaeta and Methanomassiliicoccales species. Three other smaller sub communities (module 1, 2 and 4) were identi ed with speci c taxonomic compositions.
Consistently with methanogen community structure and composition, the abundance of sub communities varied according to sampling location in Lake Remoray. Two speci c methanogen sub communities were found exclusively in the sediment of the deeper area (module 4) and sediments close to tributaries (module 1), thus further con rming the existence of diverse methanogen niches across Lake Remoray sediment and again highlighting the need to integrate lake-scale characterization of methanogen communities when assessing lake methane emissions.
Linking spatial variation of methanogen communities to methane production rate in Lake Remoray Analysis of module structure in the co-occurrence network can also provide substantial information on the functional and ecological properties of various ecosystems [90]. For instance, numerous studies have used module structure to explain the in uences of environmental parameters on microbial community structure in lake ecosystems [92][93][94][95][96]. Functional relationships between modules showed that the amount of OM in the sediments was a key factor of variation in abundance of methanogenic sub communities, in uencing four modules of the network. The quantity and quality of sediment OM are now well recognized as signi cant drivers of variation in sediment methanogen communities, shaping abundance, genetic structure and taxonomic composition [53,82,97]. In Lake Remoray, higher sediment OM content positively in uenced modules 1 and 2, which are more abundant in sediments close to tributaries. As described in section 4.1, sediment OM content and C:N ratio are higher in this part of lake, mainly due to allochthonous OM deposition by tributaries. Conversely, modules 3 and 5 were negatively in uenced by sediment OM content and their large abundances were associated to autochthonous OM depositions characterized by lower C:N ratio and related to the lake's high primary production. Taken together, these results suggest that the origin (allochthonous vs. autochthonous) and quantity of OM deposition in Lake Remoray sediment has led to a diversi cation of ecological niches for the methanogenic community.
The two larger modules (modules 3 and 5) were positively associated to methanogenesis potential and sediment methane concentration, suggesting that the methane production of Lake Remoray is mainly carried out by the biotic interactions between methanogen taxa composing these two sub communities.
However, only one of these modules was positively associated to methane production rate. On the other hand, the smaller methanogen sub community in the network (module 3) represented by the interactions between three methanogen taxa (Methanolinea, Methanofastidiosales and Methanospirillum) was strongly associated to methane production rate. This suggests that in the deeper part of Lake Remoray, methane production intensity depends mainly on the activity of a handful of methanogen species.
However, this interpretation warrants caution, as production rates were estimated as a methane potential measured under controlled in-lab conditions. The absence of any strong relationship between modules 3 and 5 and methane production rate could be due to the experimental approach used. The relationship between variation in abundance of methanogenic communities and methane production rate in lake sediment is di cult to explore [97]. Indeed, methane production rate does not rely solely on the presence and abundance of methanogenic species, and several environmental parameters also have to be factored in, such as OM composition, temperature, ion content (i.e. SO 4 2-, Fe 3+ , NO 2 -, NO 3 -) and pH [5,8,13,97].

Conclusion
Methane production and methanogenesis potential in sediments were expected to be higher in profundal areas of the lake where annual hypoxia occurs during the late summer-winter period, but we found that some littoral zones contribute as much or more to methane diffusion along the water column. Spatial patterns of methane concentration and production potential were associated with variations in nature and quality of sediment organic matter and the presence of speci c methanogenic sub communities.
Taken together, these results emphasize the importance of integrating the whole lake scale (i.e. both sediment and water column, from littoral zone to the deepest parts, including riverine zones) in order to better gauge the contribution of lake ecosystems to global methane emissions. Temporal scale is also an important factor, especially for lakes that exhibit seasonal anoxic conditions in bottom water. Further studies are needed to address the sediment-to-water column continuum and improve our predictions of methane emission from lakes. Combining molecular and geochemistry approaches with powerful isotopic approaches could valuably integrate both the spatial and temporal scales.

Funding
The authors did not receive support from any organization for the submitted work.

Con icts of interest
The authors have no con icts of interests to declare that are relevant to the content of this article.

Ethics approval
No approval of research ethics committees was required to accomplish the goals of this study because experimental work was conducted with sediment microbial species.

Availability of data and material
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Code availability
Codes for statistical analysis used in this study are available from the corresponding author on reasonable request.    Average values for methane production potential (a), methane concentration (b) and methane diffusion (c) along the rst 10 cm of sediment cores for littoral-zone (light blue) and profundal-zone (blue) areas of Lake Remoray. Symbols represent the different sampling points across Lake Remoray (see Fig. 1).
Methane diffusion (c) represents the sum of methane uxes calculated between each stratum for each sampling point. Microbial density (a), genetic potential of methanogens (b) and methanotrophs (b) along the sediment cores for littoral zone (light blue) and profundal-zone (blue) area of Lake Remoray. Histograms and vertical lines represent the averages of sampling points from littoral and profundal zones for each sampling stratum of sediment cores, with standard deviation (n=5). Symbols represent the different sampling points across Lake Remoray (see Fig. 1). Asterisks indicate signi cant differences between littoral-zone and profundal-zone sediment. Non-metric multidimensional scaling (NMDS) ordination plot derived from weighted pairwise Unifrac distances for (a) bacterial and archaeal and (b) methanogen communities for all strata of sediment cores sampled across Lake Remoray. Symbols represent sampling points, and numbers inside represent the sediment stratum. Stress values for the two ordination plots were < 0.2, which indicates that these data were well represented by the two-dimensional scaling. Vectors in the bi-plot overlay represent signi cant correlations (P<0.05) between phylum abundance and genetic structure. The angle and length of the vector indicate the direction and strength of the variable.

Figure 6
Phylogenetic tree of 41 methanogen OTUs extracted from the 16S rRNA gene sequencing dataset (a) and log2 fold changes in relative abundance of methanogen OTUs between sediment cores sampled in littoral-zone and profundal-zone areas of Lake Remoray (b). Each symbol represents the comparison between littoral-zone and profundal-zone areas for a given stratum and OTU. Numbers in white represent the stratum depths (in cm) where signi cant change was observed. OTUs not classi ed at genus level (NC) were assigned to order level.

Figure 7
Co-occurrence networks of methanogen communities in the rst 10 cm of Lake Remoray sediment cores.