Introduction

Methane is responsible for about 20% of the warming induced by long-lived greenhouse gases since preindustrial times [1]. Freshwater lakes significantly add to the natural global methane emissions [2, 3]. On the other hand methane-oxidizing bacteria (MOB) greatly reduce methane emissions from lakes to the atmosphere [4]. Typically, many MOB taxa co-occur in lake water columns [5], although all rely on the same primary resources, methane and oxygen, for growth. This raises the question of how diversity within this functional group is maintained or how the “competitive exclusion principle” [6] is bypassed. Niche partitioning could be one possible explanation for the observed diversity [7], and might further be a mechanism that increases the efficiency [8] of the methane filter in lakes.

Lacustrine MOB usually belong to Alpha- or Gammaproteobacteria, which oxidize methane with oxygen [9]. However, high abundance of Candidatus Methylomirabilis limnetica (NC10 phylum) has recently been reported in a lake water column [10]. Candidatus Methylomirabilis limnetica likely uses an intra-aerobic pathway for methane oxidation, forming oxygen from nitrite like Ca. Methylomirabilis oxyfera [11]. Other anaerobic methane oxidizers, i.e., ANME or related archaea [12,13,14], were found in lake sediments [15, 16], but seem to be of minor importance in anoxic water columns [17, 18]. Although the coarse phylogenetic affiliation of lacustrine MOB is known, several uncultivated lake and aquatic clusters within type Ia and Ib MOB (Gammaproteobacteria) are defined solely based on environmental pmoA sequences [9, 19]. Sequencing of pmoA is widely used for MOB detection in the environment [20]. However, the ecology and relevance of many uncultivated MOB remain largely unknown.

The oxygen–methane counter gradient of stratified lakes potentially provides distinct niches for MOB: methane-deficient conditions with excess oxygen, oxygen-deficient conditions with methane in excess, and a narrow interface with possibly non-limiting resource supply (Fig. 1). MOB diversity and efficient methane filtering from stratified lakes may depend on niche adaptations of MOB to ambient oxygen and methane concentrations. Previous studies indeed indicated vertical shifts in lacustrine MOB assemblages below genus level resolution [5, 21, 22]. Although the oxygen–methane interface is generally considered as the main site of methane oxidation, MOB have been shown to be present and active in oxygen-deficient water layers [23]. Fermentation and denitrification serve as oxygen-saving strategies in some isolates [24,25,26]. Further, oxygen availability has been shown to impact MOB species assemblage in incubations [27]. Isolates are known to have different affinities for methane [28], and in soils different subsets of MOB became active, depending on the supplied methane concentration [29]. In lakes, the oxygen–methane counter gradient is sometimes accompanied by a temperature gradient which might also act on the MOB assemblage [30]. In addition, other resources like nitrogen or grazing might play a role in structuring MOB niches [31]. Despite the evidence for MOB adaptations to resource availability, niche partitioning has not been demonstrated for lacustrine MOB.

Fig. 1
figure 1

Schematic placement of potential MOB niches along the vertical oxygen–methane counter gradient. a Schematic resource gradients of oxygen and methane in proximity to the oxygen–methane interface. b Potential main niches based on resource availability. c Schematic expected distribution patterns of MOB ASVs occupying differential niches

Full size image

In this study we investigated to what extent niche partitioning of MOB occurs along the natural oxygen–methane counter gradient in four stratified lakes in order to better understand MOB diversity and the structure of the effective methane filter in lakes. We focused on the section of the water column with oxygen and methane gradients above and below the oxygen–methane interface in these lakes to systematically resolve changes in the MOB assemblage and their potential activity. We compared MOB distribution patterns based on 16S rRNA and rRNA gene, as well as pmoA DNA and mRNA sequencing. By correlating MOB distribution patterns to oxygen, methane, and temperature, we were able to assign tentative niches to several lacustrine MOB at ASV (amplicon sequence variant) level, including uncultivated MOB.

Materials and methods

Sample collection and physicochemical lake profiling

We selected four eutrophic Swiss lakes known to accumulate methane in the hypolimnion. Two representatives of two common stratified lake types, seasonally and permanently stratified, were chosen. Seasonally stratified Rotsee [32] (16 m max. depth) and Greifensee [33] (18 m max. depth) were sampled during stratification on July 27th and August 20th 2015, respectively. Permanently stratified Lake Zug (198 m max. depth) and Lake Lugano (288 m max. depth) [34] were sampled October 12th and November 4th 2015, respectively. Temperature, oxygen, and pressure were profiled with a custom-made in situ analyzer (PIA) [34], equipped with a remote-controlled 50 mL syringe rosette sampler. PIA enabled sampling of steep gradients, if necessary. Compared with commonly used Niskin bottles (depth integration >50 cm), PIA integrates <5 cm allowing a vertical sampling resolution of 25 cm. Steepness and depth extension of the oxygen–methane gradient differ between lakes. To cover similar concentration ranges, vertical spacing of samples was adapted for each lake. Sampling depths were chosen on-site based on preliminary oxygen profiles from two micro-optodes (PSt1 and TOS7, Presens, Regensburg, Germany), the latter with a detection limit down to about 20 nmol L−1 [34].

Water samples for nitrate, nitrite, ammonium, sulfate, and dissolved inorganic carbon (DIC) analyses were pre-filtered with a cellulose acetate filter (0.2 µm pore size, Sartorius, Göttingen, Germany) and stored at 4 °C until analysis. Nitrate and nitrite were quantified with flow-injection analysis (ISO 13395, SAN++, Skalar, Breda, The Netherlands), sulfate with ion chromatography (881 IC compact plus, Metrohm, Herisau, Switzerland) and DIC and nonpurgeable organic carbon with a total organic carbon analyzer (TOC-L, Shimadzu, Kyoto, Japan). Nitrite was measured in Rotsee, Greifensee, and Lake Zug. Ammonium was measured photometrically according to [35]. For a summary of environmental parameters, see Supplementary Table 1.

For methane measurements, 20 ml lake water was filled into sealed 40 ml serum bottles and methane concentrations were derived from headspace concentrations measured by gas chromatography. Details are provided in the Supplementary methods.

Potential methane oxidation rates

We used a modified protocol [17, 36] to measure potential methane oxidation rates by incubations with 13C–CH4 and oxygen. Rates were calculated from the linear increase of 13C–CO2 over five time points (0, 3, 9, 21, 44 h). Details are provided in the Supplementary methods.

Nucleic acid preparation

Immediately after sampling, we filtered 110–220 ml of lake water onto 0.2 µm pore size cellulose acetate filters (Sartorius) and froze them on dry ice. Filters were stored at −80 °C until extraction. DNA and RNA were extracted simultaneously with the AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany) with bead beating (glass beads: 150–212 µm) on a FastPrep-24 homogenizer (45 sec, step 5, MP Biomedicals, Santa Ana, CA, USA). Residual DNA from the RNA was removed with the Ambion TURBO DNA-free Kit (Thermo Fisher Scientific, Waltham, MA, USA), which was checked by 16S rRNA gene PCR (27f/1492r, 35 cycles) and gel electrophoresis. If necessary, we repeated the DNA removal step. RNA was reverse-transcribed with random hexamers into cDNA with the SuperScript IV First-Strand Synthesis System (Invitrogen, Thermo Fisher Scientific), which was further used for 16S rRNA and pmoA mRNA library preparation.

Library preparation and sequencing

Sequencing of 16S rRNA gene was used to identify MOB within their well-described phylogenetic clades and to estimate their relative abundance within the bacterial community. Reverse-transcribed rRNA was used to verify protein synthesis potential. The commonly used functional marker gene pmoA [9], encoded by most MOB provided an independent assessment of the assemblage and its mRNA provided further information on MOB activity. The 16S rRNA gene/rRNA and pmoA DNA/mRNA amplicon libraries were sequenced in two runs on an Illumina MiSeq platform (Illumina Inc., San Diego, CA, USA) with 600-cycle MiSeq reagent kit v3 (Illumina Inc.) at the Genetic Diversity Centre (ETH Zurich) to obtain >10,000 reads per sample. PCR conditions and library preparation are detailed in the Supplementary methods.

Sequence analyses

We used the DADA2 workflow (version 1.6.0) [37] in R (3.4.2) [38] to determine amplicon sequence variants (ASVs) of pmoA genes and their mRNA transcripts, and 16S rRNA genes and 16S rRNA [39]. pmoA ASVs were translated to amino acid sequences (aaASVs) in MEGA7 [40], and aaASVs reaching >2% relative abundance in at least one sample of the respective lake were retained. 16S rRNA ASVs were classified based on the SILVA database (v132) [41]. ASVs affiliated with known methanotrophic groups (Order Methylococcales, Genus Methylocystis, and Genus Candidatus Methylomirabilis) were used for further analyses if reaching >0.2% relative abundance in at least two samples of the respective lake. Phylogenetic trees of MOB 16S rRNA gene and pmoA mRNA were calculated with MEGA7 neighbor-joining method using Jukes-Cantor evolutionary distance or Poisson correction method, respectively, and 10,000 bootstrap replications. Details on the amplicon sequence analysis are provided in the Supplementary methods.

Sequence accession numbers

Raw reads were submitted to the European Nucleotide Archive under project number PRJEB27579.

Quantitative real-time PCR

Copy numbers of pmoA genes and transcripts were quantified with qPCR in a 10 µl volume using SYBR Green I Master chemistry (Roche Diagnostics, Rotkreuz, Switzerland) with 0.2 µM primer (189f, mb661) [42] and 2 µl of 1:10 diluted DNA or cDNA on a LightCycler 480 (Roche). The cycling program of [43] was adapted by using 10 min initial denaturation, 15 s denaturation during cycling and data acquisition at 79 °C. The standard, a plasmid containing the pmoA gene of Methylococcus capsulatus was measured in quadruplicates at eight different dilutions (5 × 107 to five copies per reaction). Samples were run in triplicate and quantified by means of the second derivative maximum analysis method of the LightCycler 480 software (version 1.5.1.62).

Statistical analyses

Depth distributions of relative abundance of ASVs and aaASVs, normalized to maximum abundance in the profile, were hierarchically clustered based on Pearson’s distance using pvclust package (2.0.0) in R (3.5.0), UPGMA agglomeration and 10,000 bootstrap replications. Clusters representing similar depth distribution were assigned based on a distance threshold individually chosen for each lake. All clusters have a reasonable bootstrapping support of AU p values ≥80% (AU, approximately unbiased p values [44]). Venn diagrams showing taxon overlap between lakes were created with the VennDiagram package (1.6.18) in R [45]. Unconstrained detrended correspondence analysis of MOB 16S rRNA gene ASVs by lake with samples colored by niche is shown in Supplementary Fig. S1. Spearman correlation p values were corrected for multiple testing [46].

Canonical correspondence analysis (CCA) was conducted to show preferences of MOB along the gradients of oxygen, methane, and temperature across all four lakes and to test the explanatory power of the variables linked to our niche hypothesis. Constrained axes were significant (p < 0.001, anova.cca, vegan [47]). Sample scores are shown in Supplementary Fig. S2. For its ability to handle linear and unimodal taxa distributions CCA was chosen. The CCAs were performed with phyloseq [48] (1.24.2) in R, based on Chi-square dissimilarity matrix of relative abundances (sample sum = 1), and centered and scaled physicochemical variables (values < LOQ set to zero). Additional to this hypothesis-driven CCA, a CCA was performed based on a subset of environmental variables (temperature, oxygen, methane, nitrate, and sulfate concentration, Supplementary Fig. S3) curated from available environmental variables by forward selection and removal of variables with a variance inflation factor above ten. Details are provided in the Supplementary methods.

Results

Oxygen–methane counter gradients and definition of potential MOB niches

In all four lakes we found vertical oxygen–methane counter gradients; oxygen decreased with depth, whereas methane increased in the oxygen-deficient part (Fig. 2a). The oxygen–methane counter gradients sampled spanned different concentration ranges, the smallest being in Lake Zug (max. 37 µmol L−1 O2 and max. 4 µmol L−1 CH4) followed by Greifensee and Lake Lugano. The largest range was found in Rotsee with up to 400 µmol L−1 oxygen and up to 228 µmol L−1 methane – the only case with a distinct oxygen–methane overlap at the interface (Fig. 2a, Rotsee). Vertical extension and steepness of the gradients varied between lakes, therefore we sampled over 4.25 m in Rotsee (25 cm resolution), 7.2 m in Greifensee (40 cm resolution), 40 m in Lake Zug (2 or 4 m resolution), and 80 m in Lake Lugano (8 m resolution). Temperature profiles are shown in Supplementary Fig. S4.

Fig. 2
figure 2

Vertical distribution patterns of representative methanotrophic 16S rRNA gene ASVs in Rotsee, Greifensee, Lake Zug, and Lake Lugano along the oxygen–methane counter gradient. a Chemical profiles. The dashed horizontal line indicates the depth of oxygen depletion to <1 µmol L−1. Proposed niches are highlighted in orange, green, gray, and blue. Values below limit of quantification are not shown for CH4, NO3−, and NO2−. Note that nitrite was not measured in Lake Lugano. Note that in Greifensee the interface was not sampled, but the position of the niche is indicated. Mind the different y-axis scales for depth (m), and x-axis scales (concentration) for the different lakes. b Proportion of Methylococcales and Ca. Methylomirabilis limnetica of bacterial 16S rRNA gene sequences. c The most abundant ASV from each cluster in d is depicted, representing the distribution of the respective cluster. On top, the maximum relative abundance of each ASV compared with all bacterial 16S rRNA gene sequences is shown. Colors are ASV-specific. d Cluster analysis based on Pearson’s distance of the depth distribution patterns of relative abundance of MOB 16S rRNA gene ASVs. Dark gray boxes depict clusters identified. The scale shows the Pearson’s distance. Raw data of Ca. Methylomirabilis limnetica relative abundances, nitrate and methane data in Lake Zug have been used for a previously published figure in the Supporting information [10]

Full size image

Potential niches were defined based on differential resource availability. A methane-deficient but excess oxygen niche (orange background, Fig. 2) was defined as the part of the gradient from the shallowest sample to the depth where either oxygen dropped below 1 µmol L−1 (dashed line) or methane concentrations started to increase relative to samples above (Rotsee). Methane generally did not exceed ~1 µmol L−1. An oxygen-deficient and excess methane niche (blue) was designated where oxygen was below 1 µmol L−1 and methane levels increased relative to samples above. The interface in between these main zones varied in its characteristics: in Rotsee oxygen and methane overlapped at concentrations >1 µmol L−1 (green). In contrast, in the other lakes we observed a methane–oxygen gap with undetectable or very low concentrations (<1 µmol L−1) of both substrates (gray). The proposed niches have internal gradients in substrate concentrations and therefore cannot be considered as internally homogeneous habitats.

Analysis of MOB diversity by 16S rRNA gene sequencing

To identify and study MOB in the context of the general bacterial community, we first used bacterial 16S rRNA gene sequencing. This approach revealed that MOB ASVs affiliated with Methylococcales (Gammaproteobacteria) were present in all investigated lake samples, whereas Ca. Methylomirabilis limnetica (NC10 phylum) was only detected in Lake Zug and Lake Lugano samples (Fig. 2b). Summing the ASVs of both clades, the maximum relative abundance of MOB reached at a specific depth was 11% in Rotsee, 15% in Greifensee, 21% in Lake Zug, and 26% in Lake Lugano. Between lakes, both unique and shared ASVs occurred (Fig. 3a).

Fig. 3
figure 3

Venn diagrams showing number of shared a MOB ASVs based on the 16S rRNA gene and b amino acid ASVs based on pmoA mRNA between the four lakes: R = Rotsee, G = Greifensee, L = Lake Lugano, and Z = Lake Zug. ASV and aaASV identities are listed in boxes, which are colored accordingly

Full size image

Type Ia Methylococcales dominated the MOB assemblage. However, the clade consisted of 30 ASVs of several distinct genera (Fig. 4a), whose relative abundance varied considerably between different lakes and depths (Fig. 2c). Only Greifensee additionally harbored one type Ib Methylococcales (ASV_1566, Fig. 4a). Alphaproteobacterial MOB sequences occasionally occurred in the raw data, but never exceeded filtering thresholds. ASV_18 with 100% sequence identity to Ca. Methylomirabilis limnetica reached up to 6 and 11% of 16S rRNA gene reads in Lake Zug and Lugano, respectively (Fig. 2b). Currently known Ca. Methylomirabilis require nitrite to oxidize methane [11, 49]. Nitrite exceeded the limit of quantification (LOQ, 0.4 µmol L−1) at three depths in Lake Zug, but we do not have measurements for Lake Lugano (Fig. 2a, Supplementary Table 1). Nitrate concentrations decreased with depth in Lake Zug and Lake Lugano, in the latter falling below LOQ (<4 µmol L−1) in the oxygen-deficient part (Fig. 2a). Hence, nitrate-reducing microorganisms might provide nitrite.

Fig. 4
figure 4

Neighbor-joining tree of methanotroph 16S rRNA gene ASVs and pmoA mRNA aaASVs retrieved from the oxygen–methane counter gradients of Rotsee, Greifensee, Lake Zug, and Lake Lugano. a Partial 16S rRNA gene sequences (423 bp) were used to calculate a tree with the Neighbor-Joining method using Jukes–Cantor evolutionary distance. b Partial pmoA amino acid sequences (155 positions) were used to calculate a tree with the Neighbor-Joining method using the Poisson correction. Nodes with bootstrap values >0.7 (10,000 replicates) are highlighted with a black circle. Accession numbers of reference sequences are to the left of their names. The scale bars (bottom) display changes per nucleotide or amino acid position

Full size image

Depth distribution and niche preference of MOB ASVs

We next analyzed the depth distribution of the identified MOB ASVs along the oxygen–methane counter gradient in each lake. Hierarchical clustering of the depth distributions (relative abundance normalized to the maximum value) based on Pearson’s distance revealed groups of ASVs differing with respect to their distribution pattern along the counter gradient and width of their distribution, both indicating niche partitioning. Distribution patterns that conform to preferences for the potential niches can be found in all lakes (Fig. 2d). For example, in Rotsee the relative abundance of groups represented by ASV_5, ASV_35, ASV_8, and ASV_4 (Fig. 2c, Rotsee) each showed a broad, roughly unimodal distribution with depth. However, the maximum and bulk of the respective distributions centered on the methane-deficient part of the counter gradient (ASV_5), the interface with the oxygen–methane overlap (ASV_35), the interface and below (ASV_8), and the oxygen-deficient part (ASV_4), respectively. Notably, ASV_166 and ASV_98 in Rotsee showed a very narrow distribution restricted to the interface (Supplementary Fig. S5). Taxa peaking in the interface were also observed in the other lakes without overlapping methane and oxygen, but these were not interface-exclusive. In contrast, several ASVs in Lake Zug did not show a clearly defined abundance maximum or trend (ASV_39, ASV_67, ASV_4, and ASV_15, Fig. 2c).

Analysis of potentially active MOB by 16S rRNA sequencing

According to the above analysis, MOB occurred throughout the counter gradient, including parts with low-resource concentrations, raising questions regarding the metabolic state of these populations. We therefore applied 16S rRNA amplicon sequencing to RNA extracted simultaneously with DNA to assess protein synthesis potential [50]. The depth distributions of MOB rRNA ASVs appeared highly similar to the rRNA gene analysis presented above (Fig. 2c, Supplementary Figs. S5, S6), reflected by a strong correlation of summed up 16S rRNA gene and rRNA MOB relative abundances in each lake (Spearman’s rho = 0.73–0.99, p < 0.001, Table 1). Likewise, significant correlations were found between rRNA gene and rRNA based relative abundances of individual MOB ASVs in each lake (Spearman’s rho = 0.88–0.95, p < 0.001). This indicates that the observed MOB maintain their protein synthesis potential and are potentially active.

Table 1 Spearman correlation coefficients of potential methane oxidation rates and abundance of pmoA mRNA/DNA and 16S rRNA/rRNA gene-based relative abundance of methanotrophs for each lake
Full size table

Linking transcriptionally active MOB and potential methane oxidation rates

Further confirmation of activity was sought from analysis of transcriptionally active MOB. For this purpose, we analyzed pmoA mRNA by subjecting cDNA obtained from RNA extracts of all samples to high-throughput pmoA amplicon sequencing. Sequence variants of pmoA mRNA translated into 30 distinct amino acid amplicon sequence variants (aaASV) that exceeded the 2% abundance threshold. Confirming the 16S rRNA gene and rRNA results, sequences affiliated with type Ia Methylococcales were dominant and diverse, with 29 aaASVs (Fig. 4b). In accordance with the 16S rRNA gene results, a type Ib aaASV occurred in Greifensee (aaASV22, Fig. 4). The primers we used for sequencing and qPCR targeted proteobacterial pmoA; therefore, the pmoA of Ca. Methylomirabilis limnetica was not detected.

Similar to 16S rRNA gene and rRNA results, depth distribution profiles of pmoA mRNA aaASVs clustered into differential patterns (Fig. 5a, b), indicating that these populations actively transcribed pmoA. The mostly broad depth distribution of MOB populations with different depth maxima, alongside clusters with narrower depth distribution was likewise confirmed (Fig. 5b). For instance, in Rotsee aaASV1 (Fig. 5b) dominated from 8–9 m depth in the methane-deficient part, similar to ASV_5 (Fig. 2c). The distribution of aaASV2 (Fig. 5b) at the oxygen–methane interface in Rotsee resembled the distribution of ASV_35 (Fig. 2c), indicating transcriptional activity at the interface.

Fig. 5
figure 5

Vertical profiles of representative pmoA mRNA aaASVs in Rotsee, Greifensee, Lake Zug, and Lake Lugano along the oxygen–methane counter gradient. a Cluster analysis based on Pearson’s distance of relative abundance patterns of the pmoA mRNA aaASVs (unique amino acid pmoA sequences), which clusters aaASVs with similar depth distribution. Dark gray boxes depict the identified clusters. b The most abundant aaASV from each cluster is depicted representing the distribution of the respective cluster. On top, the maximum relative abundance of each aaASV compared with all pmoA mRNA sequences is shown. Colors are aaASV specific and highlight the most abundant aaASV of each cluster in a. c Corresponding profiles of pmoA mRNA copies L−1 determined by qPCR and potential methane oxidation rates (MOX) derived from 13C–CH4 incubations are shown. The dashed line indicates oxygen depletion to <1 µmol L−1. Tentative main niches are highlighted in orange, green, gray, and blue (see Fig. 1 and “Results”). Mind the different scales of the axes between lakes

Full size image

To support the link between pmoA transcript abundance and methane oxidation activity we compared them in each lake (Fig. 5c). All samples, except for the shallowest measured depth in Lake Lugano, contained pmoA transcripts as detectable by qPCR. Potential methane oxidation rates were significantly correlated with pmoA mRNA copies L−1 in all lakes (p < 0.05), except Lake Lugano (Table 1). Similarly, MOB 16S rRNA gene relative abundance was significantly correlated with potential methane oxidation rates in all lakes (p < 0.001), except Lake Lugano (Table 1). Taken together the data indicate that the observed depth distributions of both 16S rRNA gene MOB ASVs and pmoA mRNA aaASVs largely reflect the potentially active MOB assemblage, with the possible exception of the situation in Lake Lugano.

To see if pmoA transcripts deviated from gene abundance, we additionally investigated pmoA diversity and copy numbers at the DNA level. Gene copies of pmoA correlated only in Rotsee and Greifensee with potential methane oxidation rates (p < 0.01, Table 1). Distribution patterns of pmoA mRNA and DNA aaASVs were however remarkably similar overall (Supplementary Figs. S7, S8), with some noteworthy exceptions in Lake Zug and Lugano. In Lake Zug, aaASV1, and aaASV16 showed a pmoA transcript maximum at 120 m coinciding with a peak in potential methane oxidation that is not reflected in the pmoA gene distribution (Fig. 5, Supplementary Figs. S7, S8). This suggests that in some cases increased transcriptional activity (pmoA mRNA) and not increased abundance (pmoA DNA) may have contributed to the observed methane oxidation activity peaks. In Lake Lugano, several aaASVs likewise showed a mismatch between DNA and mRNA (e.g., aaASV13, aaASV1, and aaASV2). Most conspicuously, aaASV2 gene increased with depth, whereas its transcripts decreased, indicating decoupling of presence and activity.

General patterns of 16S rRNA gene and pmoA mRNA-based MOB assemblage across lakes

Canonical correspondence analyses (CCA) were conducted to show preferences of MOB along the gradients of oxygen, methane, and temperature across all four lakes and to test the explanatory power of the variables linked to our niche hypothesis, therefore other factors were not included. The CCAs explained 47% and 39% (three significant axes) of the variance in 16S rRNA gene and pmoA mRNA data, respectively, revealing general preferences of ASVs and aaASVs along the oxygen–methane gradient in agreement with our hypothesis (Fig. 6).

Fig. 6
figure 6

Canonical correspondence analysis (CCA) calculated based on a Chi-square dissimilarity matrix using relative abundances of all MOB sequences in a sample (sample sum = 1) and physicochemical parameters (T = temperature). The CCA was calculated for a 16S rRNA gene MOB ASVs and b pmoA mRNA aaASVs in 65 samples from the oxygen–methane counter gradient of the four different lakes. Colors are ASV/aaASV specific and correspond to the colors used in previous figures. The dot size visualizes the number of lakes in which the ASV/aaASV occurred

Full size image

In the following, we focus on taxa occurring in more than one lake. The CCA of 16S rRNA gene ASVs revealed that ASV_5 of the uncultivated CABC2E06 clade correlated most positively to increased oxygen concentration and temperature (Fig. 6a). Following the oxygen gradient, ASV_98 and ASV_166 correlated with oxygen, though to a lesser extent, which is in line with the observation that these ASVs occurred at or near the oxygen–methane interface. The majority of MOB ASVs clustered together in the oxygen-deficient, low temperature part of the biplot, indicative of a preference for the oxygen-deficient part of the counter gradient. This group included Ca. Methylomirabilis limnetica (ASV_18) and ASV_203, sequence identical to lacustrine Crenothrix. Both ASVs occurred in Lake Zug and Lugano, which contained low to moderate methane concentrations of up to 4 and 39 µmol L−1, respectively (Fig. 2a). Although ASV_8, related to lacustrine Crenothrix, was mainly present in the oxygen-deficient part, it correlated with methane.

The CCA of pmoA mRNA aaASVs confirmed that few MOB taxa correlated with oxygen (Fig. 6b). One of them, aaASV1 of the uncultivated Lake cluster 1, occurred in all lakes, like ASV_5 (Fig. 6a), and might therefore originate from the same organism. Five aaASVs found exclusively in Greifensee and aaASV6 strongly correlated with oxygen and temperature (Fig. 6b). Similar to ASV_166, aaASV16 appeared at intermediate oxygen levels (Fig. 6b). In line with the 16S rRNA gene CCA, the pmoA mRNA CCA showed that the majority of aaASVs preferentially occurred in oxygen-deficient, low-temperature samples (Fig. 6b). One of them, aaASV14, is sequence identical to pmoA of lacustrine Crenothrix and therefore probably originated from the same organism as 16S rRNA gene ASV_203, indicating transcriptional activity of lacustrine Crenothrix where it occurred. Another sequence related to lacustrine Crenothrix, aaASV3, correlated with methane, similar to 16S rRNA gene ASV_8 (Fig. 6). Overall, the CCA of pmoA mRNA aaASVs strengthened the evidence for ecological differentiation of MOB within the oxygen–methane counter gradient.

An additional parsimonious CCA model (Supplementary Fig. S3) incorporating further environmental variables suggested that nitrate and sulfate can help explain additional variance in the datasets (total explained variance 65% and 63% 16S rRNA gene and pmoA mRNA data, respectively) on five significant axes. Whereas the variance explained on the first two axes increased only slightly (2.9% and 6%, respectively) and did not reveal fundamentally different aspects of MOB assemblage structure.

Discussion

Our findings suggest that a structured assemblage of different MOB taxa with different ecological niche preferences collectively and effectively performs methane oxidation in stratified lakes (Figs. 2, 5). The four lakes harbored both unique and common MOB taxa along the oxygen–methane resource gradient and we discern general patterns of niche preference for some of the common MOB taxa (Fig. 6). MOB taxon distributions aligned well with oxygen and methane concentrations, indicating specialization for distinct ecological niches along the counter gradient (Fig. 2c).

Main MOB clades of the oxygen–methane counter gradients

In line with previous studies on stratified lakes [21, 51, 52], type Ia MOB dominated the MOB assemblage in terms of diversity and abundance. Different vertical distribution patterns of type Ia MOB taxa, evident from rRNA gene, rRNA, and mRNA analysis (Figs. 2c, 5b), indicated niche differentiation at a low taxonomic level as suggested in earlier DGGE and cloning-based studies on stratified lakes [5, 21, 22]. In Lake Zug and Lugano another MOB clade, denitrifying Ca. Methylomirabilis (NC10 phylum) [11], showed almost equal abundance (Fig. 2b), but consisted of only one taxon (ASV_18). Recently, the high abundance and transcriptional activity of this novel candidate species, Ca. Methylomirabilis limnetica, have been reported from Lake Zug [10]. We additionally observed Ca. Methylomirabilis limnetica in Lake Lugano, suggesting a more important role of these denitrifying MOB in permanently stratified lakes than previously thought.

The three tentative niches correspond to distinct regions of oxygen and methane availability (Fig. 1). Oxygen and methane concentrations together with temperature explained up to 47% (CCA of ASVs) of the observed variance. Although other often co-correlated parameters may act as potential co-drivers (C Guggenheim, unpublished), methane and oxygen are obvious drivers constituting the main growth substrates, which have previously been shown to act selectively on MOB taxa [29, 53]. In the following, we discuss differential distribution patterns and tentative niches of MOB taxa in the oxygen–methane counter gradient from top to bottom with a focus on taxa occurring in more than one lake.

Methane-deficient niche

The methane-deficient niche with increased oxygen concentrations was occupied by ASV_5 of the uncultivated CABC2E06 clade (Fig. 2c) and aaASV1 (pmoA mRNA) of uncultivated Lake cluster 1 [9], which might stem from the same organism. Although this CABC2E06 sequence has been detected in several aquatic environments [54, 55], its ecology and physiology has so far remained elusive. Our observations allow some hypotheses regarding its ecology to be formulated: to thrive at increased oxygen concentrations requires tolerance for elevated oxygen concentrations as oxygen sensitivity has been shown for some MOB strains [56, 57]. Methane half saturation constants of low-affinity MOB are usually >1 µM [28] and methane concentrations in the methane-deficient part of the water column are often lower than that, thus a high methane affinity would be advantageous in this habitat. High-affinity variants of particulate methane monooxygenase occur in the alphaproteobacterial Methylocystis/Methylosinus group [58, 59], but also representatives of the upland soil cluster (Gammaproteobacteria) have been linked to high-affinity methanotrophy in soils [60]. We speculate that the organism behind the ASV_5/aaASV1 pair has a relatively high methane affinity and tolerates high oxygen concentrations. Confirmation of these traits would, however, require isolation or enrichment. These traits might be shared by other ASVs e.g., ASV_1566 and aaASV22 based on their distribution pattern (Fig. 6).

Oxygen–methane interface

The oxygen–methane interface is a narrow zone with potential non-limiting concentrations of both substrates (Rotsee, Fig. 2a). MOB with fast growth rates at non-limiting methane and oxygen concentrations could thrive at the interface. The distribution patterns of ASV_166 (related to Methylomonas) and aaASV16 (related to Methylovulum), both present in Rotsee and Lake Zug, are congruent with such a niche (Supplementary Fig. S5, Fig. 5b). Intuitively, one would expect most aerobic MOB to thrive at the oxygen–methane interface, but a pronounced peak of MOB relative abundance at this depth was not observed in the lakes studied (Fig. 2b). However, at least in Rotsee and Lake Zug the aerobic MOB assemblage showed maximum potential methane oxidation rates and pmoA mRNA copy numbers close to the interface (Fig. 5c). Only Rotsee showed a pronounced oxygen–methane overlap at the interface, and only here did we observe ASVs that were strictly confined to the interface. Given this spatial narrowness and lack of a pronounced MOB abundance maximum in this zone in some lakes, it appears likely that much of the methane oxidation takes place above and below the oxygen–methane interface, proceeding under resource-limited conditions. An oxygen–methane overlap might not always be present in stratified lakes [10] due to e.g., significant anaerobic methane oxidation [10] or cryptic oxygen cycling [17, 18, 61] preventing methane from reaching the oxygen gradient. Therefore, this niche might not be present in all lakes or might be temporally restricted. Further, the oxygen–methane interface might be missed by conventional sampling, as it spanned <1 m in case of Rotsee and Greifensee. In the deep lakes where our samples were spaced more widely, we cannot exclude additional fine structure in MOB populations around the interface. However, the gradual depth distributions of many observed taxa did not indicate that we missed important fine structure.

Oxygen-deficient niche

Below the oxygen–methane interface oxygen-deficient conditions with excess methane concentrations appear. Intermittent nanomolar oxygen intrusions can occur, whose depth penetration fluctuates up to several meters within 1 day [34]. Some MOB are adapted to oxygen deficiency. Two potential pathways, fermentation and denitrification, have been shown to conserve energy in aerobic MOB isolates [24,25,26]. Ca. Methylomirabilis oxyfera-like bacteria, on the other hand, are specifically adapted to anaerobic conditions by generating oxygen internally [11]. Lacustrine Crenothrix (ASV_203/aaASV14) preferred the oxygen-deficient part of Lake Zug and Lugano (Figs. 2, 5). Lacustrine Crenothrix has previously been shown to oxidize methane under oxygen-deficient conditions and possesses genes for respiratory nitrate reductase and nitrite reductase, which might allow for survival under oxygen-deficient conditions [52]. A closely related sequence, ASV_8 (97% identity to lacustrine Crenothrix), which was present in all lakes investigated and abundant in Rotsee and Greifensee (Supplementary Fig. S5), strongly correlated with increased methane concentration (Fig. 6). Similarly, aaASV3 correlated with methane and was abundant, confirming transcriptional activity in the oxygen-deficient – high methane part (Fig. 6). In line with previous studies [23, 62], MOB affiliated with lineages that are thought to be aerobic abundantly inhabited the oxygen-deficient parts, even in the deep lakes (Lake Zug and Lugano), which remains a puzzling phenomenon [63].

Ca. Methylomirabilis limnetica was most abundant at moderate methane concentrations (4 µmol L−1 in Lake Lugano and Zug) under oxygen-deficient conditions (Fig. 2c). Ca. Methylomirabilis limnetica was absent in Greifensee and Rotsee, which might be due to photosynthetic oxygen production as light penetrates the seemingly anoxic part, likely fueling aerobic rather than anaerobic methane oxidation [17, 18]. Further, seasonal mixing in Rotsee and Greifensee might limit growth of Ca. Methylomirabilis limnetica to months with oxygen-deficient conditions (July–Nov.) [32], which might be too short to compete successfully if doubling times are in the longer range of reported values for Ca. Methylomirabilis oxyfera (1–2 weeks) [11]. In Greifensee another explanation could be the artificial aeration, which adds oxygen to the hypolimnion [64]. Although Greifensee was oxygen-deficient below 7 m depth according to our measurement, this oxygen source might contribute to the increasing MOB abundance and potential methane oxidation activity with depth (Figs. 2a, 5c).

Niche overlap

Although we found compelling evidence for niche partitioning, it is nevertheless evident that the ASVs and aaASVs substantially overlapped along the counter gradient, confirming previous findings [5]. Several ASVs clearly shared each of the described niches; some ASVs show strikingly similar depth distributions (Fig. 2d). This raises the question as to which mechanisms allow these organisms to co-occur. As shown previously for Rotsee, Lake Zug, and Lugano, the submicromolar oxygen zone fluctuates in depth, with the highest fluctuations in Lake Zug, spanning between 0.6 and 5.5 m within 1 day [34]. These fluctuations are faster than estimated environmental MOB growth rates [17], and therefore might allow for species co-occurrence and might blur niche boundaries. To what extent methane concentrations fluctuate is currently unknown as no fast-response sensors are available. An alternative explanation for co-occurrence without competition within the oxygen-deficient part could be anaerobic dormancy [65]. However, in Rotsee, Greifensee, and Lake Zug relative MOB abundance, potential methane oxidation rates, and pmoA mRNA copies were all significantly correlated (Table 1), rather contradicting the dormancy hypothesis. In Lake Lugano, methane oxidation rates did not significantly correlate with MOB relative abundance and pmoA mRNA copies, indicating that MOB dormancy may play a role in this system. Alternatively, the lower number of samples in Lake Lugano or other, so far unknown factors may explain the nonsignificant correlations.

The niche hypothesis based on methane, oxygen, and temperature certainly does not explain the entire complexity of our dataset. Some ASVs seem to change niches. For example, ASV_32, occurred in the oxygen-deficient niche in Rotsee and Greifensee, but had a substantial population in the methane-deficient niche and the interface in Lake Zug and Lugano. This could be due to microdiversity below ASV level (ecotypes) or it could indicate ecological versatility or dependence on other environmental factors not observed in our study. Moreover, we may not always cover the appropriate gradient for each ASV if their main niche was located above or below the investigated counter gradient.

In this study, we have advanced the understanding of distribution, niche partitioning, and the activity of MOB along the oxygen–methane counter gradient of stratified lakes by means of 16S rRNA gene and rRNA, and pmoA mRNA analysis, potential methane oxidation rates, and physicochemical measurements at a high vertical profiling and taxonomic resolution. We developed testable hypotheses on niche adaptations of MOB at the ASV level, including a tentative niche for a member of the uncultivated CABC2E06 clade, whose ecology is so far undefined. Future research using experimental approaches or comparative genomics could yield further insight into potential niche adaptations. We have shown that active MOB populate the oxygen–methane counter gradient and that the MOB assemblage changed from the methane-deficient to the oxygen-deficient part. This leads to the conclusion that the effective methane filter in lakes is structured, likely based on niche-specific adaptations. In addition, we propose that this niche partitioning of MOB could be a mechanism that increases their overall resource exploitation along the oxygen–methane counter gradient, as suggested for other organisms [8]. It is likely that the trait diversity of MOB makes the methane filter in lakes more effective and thus contributes importantly to reducing methane emissions from lakes.