This study is part of the larger research effort “Lake Eutrophication Impacts on Carbon Accumulations in Sediments” (LEICAS), in which the long-term consequences of eutrophication on sediment biogeochemistry, microbiology, and ecosystem ecology are investigated. LEICAS was initiated in 2016, when Lake Lucerne, Lake Zurich, Lake Zug, Lake Baldegg, and Lake Greifen were each sampled at three stations that covered a water depth gradient from shallow sublittoral to profundal sediments. Biogeochemical, organic geochemical, and 16S rRNA gene-based microbiological data on all stations were published previously by Fiskal et al. (2019, 2021a) and Han et al. (2020). An experimental study to investigate impacts of lake “oligotrophication” on sedimentary processes has also been published (Fiskal et al. 2021 b). Depth profiles of respiration zones are summarized in Table 1.

The five lakes all experience seasonal changes in primary production, thermal stratification, and sedimentation (Teranes et al., 1999; Bürgi, 2000; Naeher et al., 2013), but differ in land use and anthropogenic eutrophication histories (Fiskal et al., 2019). Despite minor increases in total P concentrations in the 1960s and 1970s, Lake Lucerne has always remained oligotrophic. Lake Zurich became eutrophic, as evidenced by the development of permanently hypoxic conditions in its deep basin, around 1890 (Naeher et al., 2013), but has been mesotrophic since ∼1980 due to improved wastewater management. The three eutrophic lakes became eutrophic at different times (Lake Baldegg: ∼1870; Lake Greifen: ∼1920; Lake Zug: ∼1930) and have remained eutrophic since. The high rates of primary production in Lake Greifen, Lake Baldegg, and Lake Zug are sustained by the retention and recycling of P that was introduced in the 20th century. Despite artificial water column mixing and aeration (Lake Baldegg: since 1982/83; Lake Greifen; since 2009), Lake Baldegg experiences strong decreases in O₂ concentrations over the summer months, and Lake Greifen continues to develop seasonal hypoxia or anoxia below 10 m water depth from June-December (for details see Supplement of Fiskal et al., 2019, and references therein).

In June and July 2016, three stations that differed in water depth were sampled in each lake (Lake Lucerne: 15, 24, 32 m; Lake Zurich: 25, 45, 137 m; Lake Zug: 25, 35, 50 m; Lake Baldegg: 21, 45, 66 m; Lake Greifen: 15, 24, 32 m) using 15-cm diameter gravity cores (UWITEC, AT). Bottom water temperatures at the time of sampling were 6–7°C at all stations except the shallow (9°C) and deep stations (5°C) in Lake Lucerne. All sites were bioturbated and contained benthic infaunal chironomids and/or oligochaetes (Fiskal et al., 2021a), except the deep, permanently hypoxic station in Lake Zurich. Sediment porewater for downstream measurements of dissolved electron acceptors (nitrate, sulfate, dissolved inorganic carbon) and respiration end products (Mn²⁺, Fe²⁺, ammonium) was sampled using rhizons (0.2 μm pore size, Rhizosphere) from a designated core with pre-drilled holes that were taped prior to coring. Solid-phase sediment samples for cell counts, DNA analyses, methane and total organic carbon (TOC) quantifications, and physical property determinations (porosity, density) were taken from an additional core using sterile cut-off syringes. Porewater cores were sampled in 1-cm depth intervals from 0–4 cm, 2-cm intervals from 4–20 cm, and 4-cm intervals from 36–40 cm. The same sampling resolution was applied during solid-phase sediment sampling, except that the top 0–2 cm were sampled at enhanced resolution (0.5 cm depth intervals). For further details on all sampling procedures, see Fiskal et al. (2019).

Additional cores were taken at each station for the measurement of O₂ concentration profiles in surface sediments using Clark-type microsensors (Unisense A/S, DK), macrofaunal community analyses, and sedimentological analyses (including Cs-137 and unsupported Pb-210 for sediment dating; further details and data in Fiskal et al. (2019, 2021a)).

O₂ and methane concentration profiles were published previously by Fiskal et al. (2019). Vertical concentration profiles of O₂ were measured using Clark-type microsensors (Unisense A/S) while methane concentrations were measured by gas chromatography using a flame-ionization detector (for details see Fiskal et al. (2019). Vertical distributions of respiration reactions, based on porewater-dissolved concentrations of O₂, nitrate, sulfate, Mn²⁺, Fe²⁺ and CH₄, were also determined in Fiskal et al. (2019).

DNA was extracted from all solid-phase samples according to Lever et al. (2015). Sediments from Lake Zug, Lake Zurich and Lake Lucerne were extracted according to lysis protocol II. Due to high humic acid concentrations that interfered with post-extraction column-based DNA purification, sediments from Lake Greifen and Lake Baldegg were treated with an additional humic acid removal step prior to purification (lysis protocol III). For detailed protocols, see Han et al. (2020).

pmoA was PCR-amplified using the general A189F - Mb661R pmoA primer combination (Holmes et al., 1995; Costello and Lidstrom, 1999) and the new NC10-pmoA239F–NC10-pmoA590R primer pair (Table 2A). The latter primer pair was specifically designed to amplify pmoA sequences of denitrifying methane-oxidizing Methylomirabilales, which are not amplified using the general primer pair but were detected in 16S rRNA gene libraries of the five lakes (Han et al., 2020). The general A189F - Mb661R combination was used both for quantitative PCR (qPCR), using complete pmoA gene sequences of Methylococcus capsulatus as qPCR standards, and for Illumina MiSeq-based community analyses.

For each of the 15 stations, qPCR analyses and sequencing were performed on 7 to 10 samples (Supplementary Table S1). Key sample targets included bottom water, sediments above and right below the oxic-anoxic transition in the upper few centimeters, samples from anoxic layers with denitrification and sulfate reduction (typically to ≤10 cm), and samples from deeper layers where methanogenesis, iron (III) and manganese (IV) reduction were the only detectable respiration reactions (Table 1). qPCR was done with SYBR Green I Master Mix on a LightCycler 480 II (Roche Molecular Systems). Each 10-μL reaction mixture consisted of 5 µl Roche Light Cycler Master Mix containing SYBR™ Green, 0.5 µl of each 50 µM primer solution, 1 µl of 1 mg ml⁻¹ bovine serum albumin, 2 µl of DNA extract, and 1 µl of molecular-grade water (for qPCR cycler protocol see Table 2B). The new NC10-pmoA239F–NC10-pmoA590R was only used for Illumina MiSeq-based community analyses but not quantification due to non-target PCR amplicons in several qPCR products.

pmoA sequencing libraries were produced according to the work flow of Han et al. (2020), which included an initial booster PCR to obtain similar amplicon concentrations across all samples, and was followed by tailed-primer PCRs (10 cycles) and index PCRs (8 cycles). For the 25-µl PCR reaction mixtures the hot start version of the TaKaRa Ex Taq DNA polymerase was used. Reaction mixtures consisted of 2.5 µl 10x Ex Taq Buffer, 2 µl dNTP solution (2.5 mM each), 0.75 µl of each primer mixture (10 µM concentrations), 2 µl of bovine serum albumin (1 mg ml⁻¹), 0.125 µl of TaKaRa Ex Taq HS, 2 µl of DNA extract, and 14.9 µl of molecular-grade water (for PCR protocols see Table 2C). Paired-end sequencing (600 cycles) was done using a MiSeq Personal Sequencer (Illumina, San Diego, CA).

All DNA sequence data were processed according to the pipeline outlined in Han et al. (2020). Briefly, raw reads were quality-checked by FastQC, followed by read end trimming by seqtk. Trimmed reads were then merged into amplicons by flash (max mismatches density, 0.15) and primer sites removed by usearch (in-silico PCR). Quality filtering was done by prinseq (GC range, 30–70; Min Q mean, 20). OTUs were then clustered according to 97% similarity thresholds (proxy for species-level) to a final number of 927 representative sequences using the USEARCH otutab (), which was also used to remove chimeric sequences. OTU sequences were then used for phylogenetic and bioinformatic analyses.

To identify potential blind spots or phylogenetic biases of the pmoA primer combinations used, we performed complementary pmoA analyses on lake sedimentary metagenomes. Metagenomes of DNA extracts from five sediment depths of the deepest station in each lake were sequenced at the Joint Genome Institute (JGI; Berkeley, CA, United States; Supplementary Table S1). Paired sequence reads were generated on an Illumina NovaSeq S4 sequencer in 2 × 150 bp mode. Assembly and annotation were performed by the standard automated pipeline of the JGI (see Supplementary Material) and loaded into the Integrated Microbial Genomes and Metagenomes (IMG/M) platform (Chen et al., 2021) for exploration and systematic gene search. Potential pmoA genes from all 25 metagenomes were identified by the AMO Pfam motif (PF02461) detected in the encoded proteins by the IMG annotation pipeline (Chen et al., 2021) resulting in 105 genes in total. The sequences were downloaded and dereplicated (100% nucleotide sequence identity) with Vsearch (v.2.14.1; Rognes et al., 2016) using standard settings, resulting in 81 unique sequences. pmoA genes were distinguished from similar amoA (ammonia monooxygenase subunit A genes) based on their placement in phylogenetic trees (see below).

Taxonomic assignments were performed using a newly constructed, comprehensive and up-to-date pmoA sequence alignment database in ARB (Ludwig et al., 2004). This database included published pmoA amplicon and metagenome sequences and was initially based on an alignment produced on a subset of sequences using ClustalW in SeaView. This gene data base has since been expanded to include >5,000 pmoA sequences in addition to representatives of pmoA-like RA21 and pxmA (putative copper-containing monooxygenase) and amoA gene clusters. In the process of database expansion, the overall alignment has been updated and optimized. All amplicon sequences obtained with the general pmoA and NC10-specific primer pairs and metagenomic pmoA sequences were added and aligned to the curated pmoA database. De-novo phylogenetic trees, calculated using ARB Neighbour-Joining with Jukes-Cantor Correction (1,000 bootstraps), were then used to taxonomically classify all pmoA OTUs to the genus-level. Tree calculations were based on a ∼470-base pair region that corresponds to the amplified region of the A189F-Mb661R primer pair.

All pmoA OTU taxonomic assignments are shown in Supplementary Material S1. Amplicon sequences are publicly accessible at the NCBI Sequence Read Archive (SRA) under Bioproject ID PRJNA785131. Samples used for metagenomes can be found under the ID Gs0142423 on JGI GOLD (https://gold.jgi.doe.gov/, Mukherjee et al., 2021). Metagenomic datasets are publicly available at the JGI genome portal (https://genome.jgi.doe.gov/portal/) under the proposal Id 504756 and at NCBI (PRJNA620348–PRJNA620356, PRJNA620360–PRJNA620364, PRJNA654670–PRJNA654680).

Diversity and richness indices were calculated with the functions “diversity” and “specnumber” in the R-package vegan (v. 2.5–7 (Oksanen et al., 2020). Kruskal–Wallis tests, combined with Dunn’s posthoc tests with Bonferroni P correction, were performed to check for significant differences in methane concentrations, pmoA copy numbers, and pmoA/16S rRNA gene ratios in relation to trophic state using R stats version 4.0.3 and the R-package FSA (Ogle et al., 2021). Non-parametric Kruskal–Wallis tests were used instead of One-way Analyses of Variance (ANOVA) due to non-normal data distributions. To test for significant relationships between MOB community composition and sediment depth, used as a proxy for geochemical gradients, and trophic state at the time of sediment deposition, Permutational Multi-variate Analyses Of Variance (PERMANOVA) were performed on genus- and OTU-level community compositions. PERMANOVAs were performed in R with the adonis2 () function in vegan based on Bray-Curtis dissimilarities (9,999 permutations). An Analysis of Similarity (ANOSIM; Bray-Curtis, 9,999 permutations) was used to test for significant differences in MOB communities between individual trophic state categories (eutrophic, mesotrophic, oligotrophic) using the anosim () function in vegan. Non-linear Multi-Dimensional Scaling (NMDS) analyses based on Bray-Curtis dissimilarities (Beals, 1984) were performed using the ‘metaMDS’ function of vegan and used to examine community dissimilarities between samples. Similarity percentage (SIMPER) analyses with the ‘simper’ function (999 permutations) in vegan were used to identify genus-level taxa and OTUs that contribute significantly to community dissimilarities between sample categories (i.e., trophic state, respiration zone). To identify potential syntrophic interactions between methanotrophic and methylotrophic microorganisms, correlations between the abundances of dominant MOB taxa and OTUs and potentially methylotrophic bacteria were determined based on Spearman’s rank correlation coefficients. The latter were calculated according to Lloréns-Rico et al. (2021) with the correlation. R script (https://github.com/raeslab/benchmark_microbiome_transformations/tree/master/scripts/, accessed on 3 July 2021).

All plots were generated with ggplot2 R-package (https://ggplot2.tidyverse.org; Wickham (2016)), basic R plotting functions, LibreOffice Calc, and Microsoft Excel.

Diffusive porewater fluxes of methane and O₂ across the sediment-water interface were calculated from vertical concentration profiles using a one-dimensional reaction-transport model (Müller et al., 2003; equation 5) applying published diffusion coefficients (D_O2 (5°C): 1.25 × 10^–5 cm² s⁻¹ (Müller et al., 2003); D_CH4 (5°C): 9.5 × 10^–6 cm² s⁻¹ (Lerman, 1979)). An empirical formation factor F of 1.2 (Maerki et al., 2004) was used to correct for porosity and tortuosity of the sediment.

Standard Gibbs energies of potential aerobic and denitrifying methanotrophic reactions, and widespread, potentially competing aerobic organotrophic and lithotrophic reactions, were calculated based on ΔG _r ^° values published in Thauer et al. (1977). The range of in situ Gibbs energies (ΔG _r ’) of these reactions in lake sediments was calculated for a range of environmentally relevant educt and product concentrations according to the equation Δ G r ' =   Δ G r ∘ + RT ln Q r where R is the universal gas constant (0.008314 kJ mol⁻¹ K⁻¹), T is standard temperature (298.15K), and Q_r is the reaction quotient. For further details, see Table 5.

Porewater-dissolved methane concentrations increased with trophic state, with sediment depth, and in several cases with station water depth (Lake Zug, Lake Zurich, Lake Lucerne; Figure 1A). Eutrophic lakes (Lakes Greifen, Baldegg, and Zug) had higher methane concentrations than the oligotrophic Lake Lucerne (p < 0.001, Kruskal–Wallis test followed by a Dunn’s post hoc test). The highly eutrophic Lake Greifen and Lake Baldegg individually also had higher methane concentrations than Lake Lucerne (p < 0.05, same test).

All stations in the eutrophic lakes and the hypoxic deep station in Lake Zurich displayed steep increases in methane concentrations from the sediment surface downward. This and the concave-down shapes of these concentration profiles indicate that net methane production sets in within the top 1 cm of sediment at these stations. These stations also had elevated methane concentrations all the way to the shallowest sampling interval (0–0.5 cm sediment depth). By contrast, methane concentrations at the middle and shallow stations in mesotrophic Lake Zurich and all stations in oligotrophic Lake Lucerne remained at lower levels (<0.5 mM) in the top cm. These stations only showed strong methane concentration increases, combined with downward concavity indicative of net methane production, below 2–4 cm.

Matching the methane concentration profiles, modeled methane fluxes increase with trophic state (Figure 2A; Table 3). The highest methane fluxes were detected at the deepest stations of all lakes except the year-round artificially aerated Lake Baldegg. By contrast dissolved O₂ fluxes into sediments did not change clearly with trophic state (Figure 2B). O₂ fluxes were consistently higher in eutrophic Lake Baldegg and Lake Zug compared to oligotrophic Lake Lucerne. The lowest O₂ fluxes were measured at the seasonally and permanently hypoxic deep stations of Lake Greifen and Lake Zurich, respectively. Ratios of methane to O₂ fluxes reflect the high methane fluxes and low O₂ fluxes at these two deep hypoxic stations (Figure 2C), where six times more (Lake Greifen) or equal amounts of methane (Lake Zurich) diffuse to the sediment surface as can be consumed by aerobic methane oxidation under the standard 1:2 stoichiometry of methane to O₂. This is not the case for the remaining 13 stations, where – assuming 1:2 stoichiometry – at most 6–30% of O₂ consumption can be accounted for by aerobic methane oxidation.

Copy numbers of pmoA, which we used as a proxy for MOB abundance, generally ranged from 10⁴ to 10⁶ g⁻¹ wet sediment and had their highest values in the top 10 cm (Figure 1B). Hereby, pmoA copy numbers varied between and within stations, but did not show consistent trends in relation to redox conditions or sediment depth. Consequently, pmoA copy numbers did not change significantly from oxic sediments to the shallowest samples below the depth of O₂ depletion (p > 0.05, here and following paragraphs: Kruskal–Wallis with Dunn’s post hoc test). Furthermore, several sites had bimodal distributions, with surface and subsurface sedimentary pmoA copy number peaks, with the latter generally found in the zone of net methane production (methanogenesis zone).

Despite having the lowest methane fluxes (Figure 2; Table 3), Lake Lucerne had the highest pmoA copy numbers (Figure 1B, p < 0.001). These were up to 10-fold higher compared to the other lakes, both within and below the O₂ penetration depth. By contrast, in Lake Baldegg, which had the highest methane fluxes (with Lake Greifen) and highest O₂ fluxes (Table 3; Figure 2), and the shallowest O₂ penetration (Table 1), pmoA copy numbers were the lowest of all lakes. These trends were statistically supported. The oligotrophic lake (Lake Lucerne) had significantly higher pmoA copy numbers than the combined eutrophic lakes (p < 0.001), but not Lake Zurich. Both Lake Zurich and Lake Lucerne also had higher pmoA copy numbers when individually compared to Lake Baldegg (p < 0.05 and p < 0.001, respectively).

The contribution of MOB to the total microbial community, expressed in % based on the ratio between pmoA and 16S rRNA gene copy numbers, ranged from ∼0.005 to 1% and was also highest in Lake Lucerne (Figure 1C). The highest values were reached at the shallow station in Lake Lucerne, despite the low methane concentrations (<2 mM; Figure 1A) and low methane fluxes (Table 3A) at this station. As for pmoA copies, Lake Lucerne had significantly higher contributions of MOB compared to the eutrophic lakes (p < 0.05), and when compared individually to Lake Baldegg (p < 0.01).

Phylogenetic analyses of pmoA sequences obtained by amplicon and metagenome sequencing revealed a phylogenetically diverse MOB assemblage (Figure 3). Groups detected include Type I and Type II MOB, Methylomirabilales, and less known pmoA-like gene clusters (RA21-like, pxma (putative copper-containing monooxygenase)).

Within the gammaproteobacterial Type I methanotrophs, diverse Type Ia and Type Ib OTUs dominated. All of these OTUs cluster more or less closely with Type Ia and Type Ib genera of the family Methylococcaceae, with the exception of several (represented by OTU_130) that cluster in the phylogenetic vicinity of Candidatus Crenothrix polyspora. Even though the latter is not phylogenetically separated from typical Type Ia Methylococcaceae pmoA (Figure 3) or 16S rRNA gene sequences (Supplementary Figure S3), it is often placed in its own family (Crenotricaceae) due to its distinct morphology. In addition to Type Ia and Type Ib pmoA sequences, Type I OTUs belonging to an unknown group (OTU_834) and the USC-gamma cluster were detected. Furthermore, one subcluster that was closely related to Methylobacter tundripaludicum was only detected by metagenome sequencing. Closest relatives of the detected Type I MOB OTUs were isolated or sequenced from diverse aquatic freshwater and marine water columns and sediments, as well as wetland and terrestrial soils.

Within the alphaproteobacterial Type II MOB, which had a lower OTU diversity than the gammaproteobacterial Type I MOB, pmoA sequences were dominated by the pmoA1 and pmoA2 clusters (collectively also known as Type IIa MOB). pmoA1 and pmoA2 encode different isoenzymes of pMMO and are found in the genera Methylocystis and Methylosinus, some members of which have both genes. In addition, we detected sequences that clustered with the genus Methylocapsa (also known as Type IIb cluster) and the less known USC-alpha cluster. Despite the known discrimination of the A1809/Mb66a primer pair against USC-alpha and Methylocapsa-like pmoA (Knief, 2015), metagenome sequencing did not yield any missing OTUs or phylogenetic clusters. The majority of Type II MOB sequences cluster with isolates and environmental sequences from terrestrial and wetland soils.

The new NC10-specific pmoA29F/590R primer combination detected sequences that cover the known phylogenetic diversity of NC10-pmoA. The fact that NC10 were also detected by metagenome sequencing, including one sequence with 100% sequence identity to Candidatus Methylomirabilis limnetica, indicates that this group occurred in significant numbers in certain samples (discussed more later). Detected NC10-pmoA sequences clustered with ones that were previously recovered from water-logged soils and aquatic sediments.

Two additional, pmoA-like clusters (RA21-like, pxma) were only detected by metagenome sequencing, matching past observations that the general A1809/Mb66a misses both of these divergent, functionally poorly understood pmoA clusters. The recovered sequences diverge strongly from previously published RA21-like or pxma metagenome sequences and are thus of unknown taxonomic origin.

The methanotrophic community detected with the general A1809/Mb66a pmoA primer pair showed a similar OTU richness across the five lakes (Figure 4A; p > 0.05, Kruskal–Wallis test followed by a Dunn’s post hoc test). The diversity and evenness was, however, significantly lower in oligotrophic Lake Lucerne compared to each of the other lakes (Figures 4B,C; p < 0.001, same tests). This indicates that, despite harboring a comparable total number of OTUs, MOB communities in Lake Lucerne are dominated by few highly abundant OTUs. By contrast, the other four lakes shared a similar, higher evenness and diversity of OTUs. Richness and diversity comparisons were not performed on NC10, as only 19 OTUs belonging to this phylum were recovered, most of which were unique to sediments of Lake Lucerne.

pmoA taxonomic compositions based on the general A1809/Mb66a pmoA primer pair showed systematic trends in relation to sediment depth and to a lesser degree trophic state (Figure 5). Clear trends within lakes in relation to water depth or from oxic to anoxic sediment layers within the same sediment cores were absent.

Type I methanotrophic communities in all five lakes were generally dominated by three taxonomic groups, termed Type Ia Unclassified Methylococcaceae, Type Ib Methylococcus, and Type Ib Unclassified Methylococcaceae (Figure 5). Bottom water was dominated by Type Ia Unclassified Methylococcaceae (in most cases >95% of reads). Exceptions are the deep stations in Lake Zug and Lake Lucerne, which are dominated by Type Ia Methylobacter and Type Ib Unclassified Methylococcaceae pmoA sequences, respectively. Communities shift systematically from bottom water into sediments, where pmoA sequences that cluster with the genus Methylococcus (Type Ib), which were nearly absent from bottom water, emerge as a second dominant group.

Percentages of Type Ia Unclassified Methylococcaceae are on average higher in sediments of today’s eu- and mesotrophic lakes compared to oligotrophic Lake Lucerne. Percentages of this group also decrease gradually with sediment depth at most stations, although these decreases do not closely match depth distributions of O₂ or anaerobic respiration reactions, or past shifts in trophic state. Similar depth- or trophic state-related changes in Type Ib Methylococcus-type pmoA contributions are less evident, though, notably, this cluster drops in relative abundances in deep sediment layers of oligotrophic Lake Lucerne and most deep oligotrophic layers of today’s eutrophic lakes. The third major Type I cluster (Type Ib Unclassified Methylococcaceae) fluctuates in percentages without clear trends in relation to O₂ or respiration zone, sediment depth, past or present trophic state. Other Type I clusters that accounted for significant, but minor percentages of pmoA reads included Type Ia Methylomonas- and Methylobacter-like sequences, Type Ib Methylocaldum-like sequences, and Unclassified Type I-sequences (Figure 5). Read percentages of these groups also did not show clear trends in relation to O₂ or respiration zone, sediment depth, past or present trophic state.

The general decreases in Type I pmoA sequence percentages with sediment depth coincided with strong increases in facultative Type II methanotrophs that mainly grouped with the Methylocystis-type pmoA1 cluster (Figure 5). This cluster was almost absent from bottom water but – independent of trophic state - became dominant below the depths of O₂, nitrate, and sulfate depletion in iron-reducing and methanogenic sediment layers.

We checked whether sedimentary trends in pmoA sequences were reproducible in the complementary 16S rRNA gene data of Han et al. (2020), and how MOB that were missed by the general A1809/Mb66a pmoA primer pair compared in relative abundances (Supplementary Figure S1). These analyses confirm the shift in dominance from gammaproteobacterial to alphaproteobacterial MOB with increasing sediment depth (and age) observed in pmoA sequence compositions. Moreover, in sediments of Lake Lucerne, members of Candidatus Methylomirabilia (NC10) dominate shallow sediments along with Type Ia Methylococcaceae. Methylomirabilia furthermore account for up to ∼30% of methanotrophic 16S reads in deep sediment layers of the eutrophic lakes, and even dominate sediments below 5 cm depth in Lake Zurich. All Methylomirabilia-related 16S rRNA gene sequences could be assigned to the order Methylomirabilales. Most reads in Lake Lucerne belong to the genus Methylomirabilis, while most reads in the other lakes group with the uncharacterized Sh765B-TzT-35 cluster.

A striking difference to the pmoA analyses was the much lower percentage of Type Ib in the 16S rRNA gene data set. In addition, alphaproteobacterial Type II were dominated by 16S rRNA sequences that cluster with the mostly methylotrophic genus Methyloceanibacter (family: Methyloligellaceae), which only includes one confirmed methanotrophic member that, however, lacks pmoA (M.oceanibacter methanicus; Vekeman et al., 2016). By comparison, percentages of sequences related to facultatively methanotrophic Methylocystis are much lower, though these also increase with sediment depth. Lastly, 16S rRNA gene data indicate the presence of low read fractions (∼0.1%) of Type III methanotrophs (Methylacidiphilaceae; Verrucomicrobia) in surface layers of the middle stations in Lake Baldegg and Lake Lucerne.

A Non-linear Multi-dimensional scaling (NMDS) ordination of samples based on major taxa recovered with the general A1809/Mb66a pmoA primer pair confirms the main inference from Figure 5 that MOB communities at the genus-level are more strongly structured in relation to sediment depth than trophic state (Figure 6A). This is confirmed based on the results of a PERMANOVA (Table 4). These show that while both sediment depth and trophic state have a significant impact on MOB community structure at the genus-level (both p < 0.001), depth explains a higher portion of the observed community variation (R ² = 0.22) than trophic state (R ² = 0.16). Complementary ANOSIM analyses, moreover, show that trophic state at the time of deposition only shows significant differences between samples that were deposited under oligotrophic vs. eutrophic conditions (R = 0.30, p < 0.001), but not between eutrophic and mesotrophic (R = 0.08, p > 0.05) or mesotrophic and oligotrophic conditions (R = -0.16, p > 0.05).

The pattern shifts when NMDS analyses are performed at the OTU-level (“species-level”; Figure 6B). Here trophic state is more strongly correlated with community structure than sediment depth, as indicated by the clear separation of oligotrophic Lake Lucerne from the other four lakes. By contrast, OTU compositions of sediment samples that were deposited under mesotrophic conditions in Lake Zurich are highly similar to samples deposited under eutrophic conditions in Lake Zurich or the other presently eutrophic lakes. Interestingly, OTU samples from the deep oligotrophic layers in Lake Zug and one deep oligotrophic layer in Lake Greifen cluster most closely with oligotrophic communities from Lake Lucerne (Figure 6B), suggesting that trophic state at the time of deposition is the key driver of MOB community structure at the OTU-level. This is confirmed based on a PERMANOVA test (Table 4), which shows that while both sediment depth and trophic state have a significant impact on MOB community structure at the OTU-level (both p < 0.001), trophic state explains a much higher portion of the observed community variation (R ² = 0.23) than sediment depth (R ² = 0.08).

Further OTU-level analyses, reveal that the three main pmoA clusters (Figure 5; Type Ia unclassified Methylococcaceae, Type Ib Methylococcus, Type II Methylocystis) in eu- and mesotrophic lakes are dominated by a larger number of OTUs than in Lake Lucerne (Supplementary Figure S2). Most of these dominant OTUs in the eu- and mesotrophic lakes are shared among these lakes. By contrast, as already indicated by Figure 4, Lake Lucerne is dominated by a smaller number of OTUs with very high relative abundances. These dominant OTUs in Lake Lucerne represent a subset of the OTUs that dominate sediments of the eu- and mesotrophic lakes. Additional similarity percentage (SIMPER) analyses reveal that some OTUs are highly abundant in oligotrophic samples of Lake Lucerne and deep oligotrophic layers of Lakes Greifen, Baldegg, and Zug, but nearly absent from sediment layers that were deposited under eutrophic conditions (e.g. Methylococcaceae OTU2, Methylocystis‐related OTU3; Supplementary Figure S2).

Abundance analyses of pmoA copy numbers and pmoA/16S ratios in relation to respiration reactions and trophic state confirm both the absence of a clear pmoA copy number relationship with respiration zone as well as the clear relationship between pmoA copy number and trophic state (Figure 7A). On average, pmoA copy numbers increase by a factor of 8 and pmoA/16S ratios increase by a factor of 4 from sediment samples that were deposited under eutrophic conditions to ones that were deposited under oligotrophic conditions. Both increases are significant (p < 0.001, Kruskal–Wallis followed by Dunn’s post hoc test).

In addition, dominant pmoA taxonomic groups (Figure 5) and their most abundant OTUs show clear patterns in relation to vertical respiration zone and to a lesser degree trophic state (Figure 7A). The observed patterns in relation to respiration zone explain the observed significant correlations of pmoA community structure at the genus- and OTU-level with sediment depth, which per se is not a likely community driver (see previous section). The Type Ia MOB genera Methylobacter and Methylomonas and Type Ib genus Methylococcus, and dominant OTUs within all three genera, have their highest calculated average pmoA abundances in oxic sediments. By contrast, the Type Ia Unclassified Methylococcaceae and Type II Methylocystis, and several dominant OTUs within these groups, have peak pmoA abundances in denitrifying, sulfate- and manganese-reducing (Unclassified Type Ia) and in manganese- and iron-reducing and methanogenic layers (Methylocystis) (Figure 7A). With respect to abundances in relation to trophic state, we generally observe that Type Ia Methylobacter and Type II pmoA sequences, and their dominant OTUs, reach their peak abundances under oligotrophic conditions, whereas Type Ia Methylomonas and Unclassified Methylococcaceae, and their dominant OTUs, reach their highest abundances under meso- and eutrophic conditions. Highly similar patterns to those in Figure 7A were observed when calculated absolute abundances of dominant taxa and OTUs were substituted with relative abundances (Supplementary Figure S4).

Calculated pmoA abundances of major MOB taxa and their dominant OTUs correlate strongly with major methylotrophic taxa and their dominant OTUs (Figure 7B). Type Ia Unclassified Methylococcaceae, Type Ia Methylomonas, and Type Ib Methylococcus, as well as dominant OTUs within these three groups, show strong positive Spearman’s correlations (Rho>0.5, p < 0.001) with 16S rRNA gene-based abundances of Methylophilaceae and dominant OTUs of Methylophilaceae. In addition, matching their observed pmoA peak abundancess in deeper layers, pmoA abundances of Methylocystis and one of its dominant OTUs (OTU_3) are correlated with 16S gene-based abundances of Methyloceanibacter and one of its dominant OTUs (OTU_45). Nonetheless, pmoA OTUs with highest abundances in the oxic or denitrifying sediment layers (Type Ia and Type Ib MOB) were much more likely to correlate with methylotrophic OTUs than pmoA OTUs reaching peak abundances in deeper anoxic layers (Type II). Among pmoA OTUs with peak abundances in the oxic or denitrifying layer, 19 and 22%, respectively, showed a significant positive correlation (Rho>0.5, p < 0.05) with methylotrophic 16S OTUs. This proportion dropped to 10 and 9% for OTUs most abundant in the sulfate and manganese reduction zones, and was 0 and 2%, respectively, in the iron reduction and methanogenesis zones. Moreover, of the 64 Methylocystis OTUs, it was only the aforementioned OTU_3 (Figure 7B) that showed significant correlations with Methyloceanibacter.

We found no clear relationship between trophic state and the percentage of pmoA OTUs that correlate significantly with the percentage of methylotrophic 16S OTUs (Rho>0.5, p < 0.05). Based on calculated absolute abundances, 10 and 12% of pmoA OTUs with peak abundances in sediments deposited under eutrophic and oligotrophic conditions, respectively, were significantly correlated with methylotrophs.

Matching measured methane concentrations (Figure 1A) and modeled methane production profiles (Fiskal et al., 2019), we observe clear increases in surface sedimentary methane fluxes in eutrophic lakes (Figure 2A; Table 3). Yet, gene copy number-based estimates indicate that MOB are on average 8-fold more abundant and account for 4-fold higher fractions of the total microbial community in samples deposited under oligotrophic compared to eutrophic conditions (Figure 7). An insufficient supply of O₂ can be ruled out as a first-order explanation as, except at the hypoxic deep stations in Lake Greifen and Lake Zurich, only a minor fraction (6–30%) of the diffusive flux of O₂ can be accounted for by aerobic methane oxidation assuming the standard 2:1 stoichiometry of O₂:CH₄ (Figure 2C).

Matching previous studies on lake sediments (Rahalkar et al., 2009; Pester et al., 2011; Tsutsumi et al., 2011; He et al., 2021; Lyautey et al., 2021), pmoA copy numbers do not change significantly from oxic to anoxic sediments (Figure 1B). Based on calculated pmoA gene copy numbers of individual MOB genera, certain Type II MOB taxa (mainly Methylocystis) even reach peak abundances in deep methanogenic and metal reducing sediment (Figure 5; Supplementary Figure S1). The observation of high or even elevated pmoA copies below the oxic-anoxic transition is at odds with cultivation data, which indicate that the vast majority of alpha and gammaproteobacterial MOB, including all Methylocystaceae, are obligate aerobes (Bowman, 2014). This raises questions concerning whether these pmoA sequences belong to catabolically active, dormant, inactive, or even dead cells.

While Type II MOB dominate deeper layers, Type I MOB (Methylococcacea) are the major group in oxic and suboxic bottom water and surface sediments (Figure 6), matching other studies on freshwater surface sediments (e.g., Rahalkar et al., 2009; Lyautey et al., 2021). Hereby the (near) absence of Type Ib, but not Type Ia, MOB in overlying water suggests a habitat preference of Type Ib for sediments. Similar to previous studies (Yang et al., 2016; He et al., 2021; Lyautey et al., 2021), Methylococcus-type pmoA dominate Type Ib communities. By contrast, Methylobacter-type Type Ia sequences, which frequently dominate lake sediments elsewhere (Tsutsumi et al., 2011; Reim et al., 2012; Yang et al., 2016), are less abundant. Instead, the majority of Type Ia sequences fall outside of known genera (‘Type Ia Unclassified Methylococcaceae’). Some of these sequences cluster loosely with Candidatus Crenothrix polyspora (Figure 3), raising questions regarding their identity as filamentous Crenotrichaceae or unicellular (mostly coccoid) Methylococcaceae. Phylogenetic trees based on both pmoA and 16S rRNA gene sequences, in which Crenotrichaceae (Crenothrix polyspora) group with well-known Methylococcaceae isolates (Figure 3, Supplementary Figure S3) could not resolve this issue. Besides Type I and Type II MOB taxa, we detect NC10-type denitrifying methanotrophs (Methylomirabilota) based on 16S rRNA gene analyses (Supplementary Figure S1). Members of the type genus Methanomirabilis occur at similar relative abundances to Type I methanotrophs in denitrifying sediments of Lake Lucerne, whereas uncultivated groups of Methylomirabilota (mainly Sh765B−TzT−35) occur at significant abundances in deeper, nitrate-depleted, fully anoxic layers of most lakes.

While the drivers of (uncultivated) Methylomirabilota communities remain uncertain, Type I and Type II MOB taxa at the genus-level show vertical distributions that match distributions of major respiration reactions (Figure 7A and Supplementary Figure S4). This raises the possibility that distributions of MOB genera in lake sediments reflect genus-level differences in the ability to live without O₂ and respire different anaerobic electron acceptors. Moreover, OTU-level analyses indicate lake trophic state as a major driver of MOB diversity and community structure at the species-level (Figure 4, 6B & 7A). Potential explanations are that differences in anabolic nutrient regimes or micro-scale habitat heterogeneity impact MOB communities at the species-level. By contrast, based on correlation analyses of MOB taxa and potential methylotrophic partner organisms (Results text), we find no evidence that co-dependencies of MOB on syntrophic partners vary with trophic state. Instead, correlations between MOB and methylotrophic bacterial OTUs follow clear taxonomic and environmental trends, being most frequently observed within Type I taxa that dominate oxic and suboxic surface sediments.

The order of magnitude higher pmoA copy numbers in oligotrophic Lake Lucerne compared to the eutrophic lakes (Figure 1B & 7A), show parallels with quantitative data on total microbial community size. Total 16S rRNA gene copy numbers were also highest in Lake Lucerne (Han et al., 2020) and were generally supported by microscopic cell counts which indicated higher microbial cell abundances in Lake Lucerne and Lake Zurich than in the three eutrophic lakes (Fiskal et al., 2019). These trends seem paradoxical given that the higher deposition and burial rates of organic carbon in the eutrophic lakes sustain higher total respiration and methanogenesis rates throughout cores from these lakes (Fiskal et al., 2019). A possible explanation is that microbial cell-specific energy requirements are higher in eutrophic lake sediment, e.g. due to enhanced physiological stress from elevated toxic contaminant loads. Yet, though more research is needed, previous studies suggest minor effects of anthropogenic toxin contamination on sediment microbial communities (Chen et al., 2014; Di Cesare et al., 2020). Notably also, the shallow and medium station in Lake Zurich are located in a Zinn (Sn)-polluted area (Lennartz, 2010), and yet had high cell counts and high qPCR values of pmoA and total 16S genes.

An alternative explanation is that – matching the higher microbial activity - sediments of eutrophic lakes have higher microbial growth rates, but that these coincide with higher rates of mortality, e.g. due to viral lysis or predation. Though we lack data on viral lysis rates, previous work has shown high abundances of oligochaetes (Tubificidae) to 20 cm sediment depth in the three eutrophic lakes and absence of these worms from all but the sediment surface in Lake Lucerne (Fiskal et al., 2021a). Gut microbiome analyses suggest protein-rich diets that are consistent with microbial grazing as a major food source of tubificids (Fiskal et al., 2021a; de Valk et al., 2017). While lake sedimentary invertebrates are known to stimulate growth of MOB in surface sediments (“microbial gardening”), leading to MOB becoming a major food source of these invertebrates in some lakes (Kiyashko et al., 2001; Deines et al., 2007; Jones and Grey, 2011), none of the stations in this study show evidence for MOB gardening (Fiskal et al., 2021a). We thus propose that non-selective grazing by tubificids causes higher mortality and ultimately results in lower abundances of microorganisms, including MOB, in sediments of the eutrophic lakes.

While non-selective grazing by oligochaetes explains lower absolute abundances of MOB, it does not explain the lower relative abundances (percentages) of MOB in the eutrophic lakes (Figure 1C, 7A, and Supplementary Figure S2). These lower percentages are particularly striking given the steep, overlapping O₂ and CH₄ concentration gradients in the top 0.5 cm of eutrophic lakes (Figure 1A), which imply high potential for aerobic methanotrophy. In the following paragraphs, we evaluate three potential hypotheses to explain the lower percentages of MOB in eutrophic lake sediments with high methane fluxes.

Direct observations on sediment cores from the eutrophic lakes revealed the presence of gas bubbles starting just below the 1 to 2 cm thick, flocculent surface layer. While gas bubbles also appeared in sediment of Lake Lucerne and the shallow and medium stations in Lake Zurich, these bubbles were restricted to deeper, more compacted sediment layers and less numerous in these locations. Based on these observations, it is possible that most methane in sediments of eutrophic lakes is present in gas bubbles rather than bioavailable dissolved form. These methane bubbles may bypass the “biological methane filter” in oxic surface sediments by rapid transit through the flocculent surface layer followed by ebullition into overlying water, explaining why methane oxidation only supports a small percentage of the microbial community in eutrophic lake sediments. The importance of ebullitive release of methane from lake sediments to overlying water and the atmosphere is well-known and especially high for shallow and/or eutrophic lakes (Thebrath et al., 1993; Schubert et al., 2012), though little is known about how ebullitive losses quantitatively compare to biological consumption by MOB in lake sediments.

Given the high sedimentation rates of labile microalgal biomass (Han et al., 2020) and shallow depletion of O₂ and nitrate in the three eutrophic lakes (Table 1), competition for O₂ or nitrate between MOB and other microorganisms could also explain the lower MOB percentages in eutrophic lakes. Thermodynamic calculations (Table 5) indicate that aerobic organotrophic microorganisms have substantially higher energy yields than aerobic or denitrifying MOB. Organotrophs could thus outcompete MOB for O₂ and nitrate, and thereby suppress methane consumption rates and the efficiency of the biological methane filter in the eutrophic lakes. By comparison, sediments of Lake Lucerne have deeper O₂ and nitrate penetration (Table 1) and receive less labile algal organic matter than the eutrophic lakes (Han et al., 2020). Here competition for O₂ or nitrate could be less intense due to stronger electron donor limitation of organotrophs, and lead to MOB being major O₂ consumers. In Lake Lucerne, MOB may even have energetic advantages over competing lithotrophic nitrifying bacteria (Table 5), which are major microbial community members in this lake (Han et al., 2020). Reduced competition for O₂ and nitrate could also explain why MOB effectively deplete upward-diffusing methane in the surface centimeters of Lake Lucerne and the oxic shallow and medium stations of Lake Zurich.

The third hypothesis concerns potential changes in dominant methanotrophic reactions in relation to trophic state. Besides performing the complete oxidation of methane to CO₂, MOB are known to under certain conditions only partially oxidize methane to low molecular weight intermediates, such as methanol, H₂, formate, or acetate, which are then utilized by denitrifying partner organisms (Kalyuzhnaya et al., 2013; Beck et al., 2013; Oshkin et al., 2014; Hernandez et al., 2015). The factors that determine whether MOB perform complete or partial oxidation are not understood. O₂ concentrations may affect the MOB genera involved, but not necessarily the occurrence of syntrophic reactions (Beck et al., 2013; Hernandez et al., 2015; Oshkin et al., 2014). Compared to the complete oxidation of methane to CO₂, the Gibbs energies of partial oxidation reactions are ∼6-fold (methanol) or 3- to 4-fold (formate, acetate) lower (Table 5). Higher contributions of energetically less favorable partial methane oxidation reactions in the eutrophic lakes could thus contribute to the lower MOB percentages in the latter. Yet, while we observe strong correlations between calculated pmoA copy numbers of certain Type I genus-level taxa and OTUs with calculated 16S gene copy numbers of Methylophilaceae (Figure 7B), the percentage of taxa or OTUs with such correlations does not change in relation to trophic state.

In conclusion, based on available data, enhanced oligochaete grazing is a plausible reason for the lower absolute abundances of MOB in eutrophic lakes, whereas the lower percentages of MOB in eutrophic lakes can be explained with enhanced methane ebullition and/or outcompetition of MOB for O₂ by aerobic organotrophs. By contrast, based on our correlation analyses, enhanced energy-sharing with syntrophic partner organisms is less likely to drive the lower percentages of MOB in eutrophic lakes.

Analyses of pmoA genus-level taxa distributions indicate strong overlaps in community structure between different trophic states, and thus argue against the existence of distinct oligotrophic or eutrophic MOB taxa (Figure 5 & 6A). Instead, the highly significant correlation of sediment depth with taxon-level MOB community structure, combined with the clear taxon-specific differences in distributions relative to vertical respiration zones, point to respiration zone as a key driver of MOB community composition at the genus-level (Figure 5, 6A, and 7A).

Interestingly, the Type Ia genera Methylobacter and Methylomonas and Type Ib genus Methylococcus, and dominant OTUs within these genera, had their highest calculated average pmoA abundances in oxic sediments (Figure 7), even though at least members of Methylobacter and Methylomonas have been linked to facultative anaerobic metabolism with or without partner organisms (Methylomonas: NO₂ ⁻, NO₃ ⁻, ferrihydrite (Kits et al., 2015b, 2015a; Zheng et al., 2020)); Methylobacter: NO₃ ⁻ (Smith et al., 2018; van Grinsven et al., 2020a)). Other MOB, including the Unclassified Type Ia Methylococcaceae have peak pmoA abundances in suboxic layers (Figure 7), suggesting that members of this paraphyletic group might couple methane oxidation to nitrate reduction. Syntrophic interactions between unclassified Type Ia Methylococcaceae and Methylomonas and methylotrophic Methylophilaceae with denitrifying potential were previously proposed (Beck et al., 2013; Cao et al., 2019; van Grinsven et al., 2020b) and would match found correlations between pmoA abundances of both groups and 16S rRNA gene-based abundances of Methylophilaceae (Figure 7B).

While the observed distributions of Type I taxa are generally consistent with previous studies, the observed trends in Type II taxa are more difficult to reconcile. Type II alphaproteobacterial Methylocystaceae sequences dominated deep anoxic sediments of all lakes (Figure 6). It remains unclear whether these deep Methylocystis-related pmoA sequences derive from metabolically active cells, resting (dormant) or inactive cells, or dead organisms. For this DNA to be from dead cells, DNA of Methylocystis would need to be more effectively preserved than DNA of Type I methanotrophs, e.g. via shielding by more degradation-resistant cell walls. Members of Methylocystis produce lipoidal cysts as resting stages, which are compositionally distinct from the carbohydrate-rich (“Azotobacter-type”) cell walls of cysts produced by most Methylococcaceae (Bowman 2014). These lipoidal cysts may be more degradation-resistant in lake sediments than Azotobacter-type cysts. Yet, enhanced stability of resting stages, which requires substantiation, does not explain why dominant Methylocystis-related OTUs reach their highest gene copy numbers in anoxic layers where respiration is dominated by metal reduction and methanogenesis (Figure 7A), or why previous studies have reported increases in relative abundances and active assimilation of methane-C by Type II methanotrophs in anoxic sediment (Yang et al., 2019; He et al., 2012). A potential explanation is that Type II Methylocystis-related MOB oxidize methane via strictly anaerobic associations with partner organisms. Yet, only one Methylocystis OTU (out of 64) showed a significant correlation with a methylotrophic bacterial OTU (genus Methyloceanibacteria; Figure 7B and related text). Thus, deep Type II MOB might depend less on syntrophic interactions than their Type I counterparts in shallower layers, or not at all. Future research will reveal whether members of Methylocystis are involved in methane oxidation, e.g. coupled to metal reduction, as was previously shown for members of Methylosinus, a sister genus of Methylocystis (Zheng et al., 2020), or even subside on substrates other than methane, such as acetate or ethanol as was previously shown for (however aerobic) members of this genus (Haque et al., 2020; Oshkin et al., 2020).

While pmoA-based community structure at the genus-level was similar across different trophic states, analyses at the OTU-level indicate trophic state as a major driver of MOB community structure at the species-level (Figure 6B). Sediments of the oligotrophic Lake Lucerne have a significantly lower OTU diversity and evenness than sediments of the meso- and eutrophic lakes (Figure 4) and are dominated by a small number of OTUs that are also present though not dominant in eu- and mesotrophic lakes (Supplementary Figure S2). Interestingly, as previously demonstrated for 16S rRNA gene-based analyses of bacterial communities (Han et al., 2020), the trophic state at the time of sediment deposition appears to better predict the present-day MOB community structure than the present-day trophic state of the lake in numerous cases. This is indicated by the fact that OTU community fingerprints of several deep samples from Lake Zug and Lake Greifen, that were deposited under oligotrophic conditions, cluster with samples from the always oligotrophic Lake Lucerne (Figure 6B). This legacy effect was most evident in a Type Ib Methylococcaceae-related (OTU_2) and a Type II Methylocystis-related OTU (OTU_3). These OTUs dominated sediment samples from Lake Lucerne and had high relative abundances in deep oligotrophic samples from Lake Zug and Lake Greifen, while being nearly absent from all other samples (Supplementary Figure S2).

Our results indicate that MOB community structure is controlled by different environmental variables on different phylogenetic levels. The strong community structuring of MOB genus-level taxa in relation to respiration zones could reflect differences in electron acceptor use or methylotrophic partners among these taxa. By contrast, the clear structuring of OTU-based communities relative to trophic state suggests that additional factors control MOB communities at the species or strain-level, such as differences in enzyme and transporter affinities in relation to nutrient regimes or competition with other microorganisms for O₂. The low community alpha-diversity and evenness in oligotrophic Lake Lucerne, where organic matter is more degraded and contains lower fractions of algal organic carbon, could reflect lower geochemical niche diversity, reduced availability of biosynthetic nutrients and reduced competition for O₂. In comparison, higher availability of nutrients needed for biomass production and higher diversity of microniches with distinct redox conditions could support higher numbers of co-existing OTUs and higher evenness in eutrophic lake sediments.