香京julia种子在线播放

Rotsee is a small (0.5 km²), eutrophic lake located in central Switzerland with a maximum depth of 16 m. Its wind-shielded location allows stable stratification from spring until mid to late autumn with an oxycline usually developed between 6 and 9 m depth. During stratification methane accumulates in the anoxic water column and reaches concentrations up to 1 mM (Schubert et al., 2010; Oswald et al., 2015). Sample collection was conducted close to the deepest point of the lake during 3 consecutive years at the beginning of stratification (June 2013), during peak stratification (August 2013), and shortly before the lake overturns (September 2014 and September 2015). Detailed methods of physico-chemical profiling, sampling, and analysis are reported in the Supplementary Material. Our physicochemical dataset includes the following parameters: conductivity (Cond), turbidity (Turb), depth (pressure), temperature (T), pH, photosynthetically active radiation (PAR, herein equated as light), concentrations of oxygen (O₂), chlorophyll a (Chl-a), total sulphide (S_Tot = H₂S, HS⁻, S²⁻), dissolved organic carbon (DOC), total dissolved nitrogen (TDN), dissolved inorganic carbon (DIC), nitrite (NO₂), nitrate (NO₃), ammonium (NH₄), sulphate (SO₄), phosphate (PO₄), dissolved (M_Diss) and total (M_Tot) metal concentrations [copper (Cu), iron (Fe), manganese (Mn), zinc (Zn), chromium (Cr)], methane (CH₄) as well as ¹³C/¹²C isotopic ratio of CH₄ (δ¹³C-CH₄). Particulate metal (M_Part) concentrations were obtained by subtracting dissolved metal from total metal concentrations. Bioavailable metal fractions (M_DGT) were measured via the Diffusive Gradients in Thin film (DGT) technique (Davison, 2016). Detailed information on the retrieval of the DGT data can be found in (Guggenheim et al., 2019).

Bacterial DNA was obtained by filtration of water samples and extraction from 0.2 μm polycarbonate filters using the PowerWater® DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA). Further details on DNA sampling and processing are reported in the Supplementary Material. The 16S rRNA and pmoA genes served as marker genes for the total bacterial community and MOB detection, respectively (Dumont, 2014). Quantitative polymerase chain reaction (qPCR) on the pmoA gene was conducted using the primer pair A189f and mb661r (5′-GGNGACTGGGACTTCTGG-3′, 5′-CCGGMGCAACGTCYTTACC-3′, Eurofins Genomics, Ebersberg, Germany), which covers most alphaproteobacterial and gammaproteobacterial MOB (Costello and Lidstrom, 1999). Although this primer set is specific for pmoA, it disfavors alphaproteobacterial MOB (Bourne et al., 2001). Verrucomicrobial pmoA sequences and sequences of NC10 phylum were not covered by the applied primer pair. pmoA amplicon libraries were prepared using the above-mentioned primers with Illumina Nextera overhang sequences at the 5′-end (Microsynth AG, Balgach, Switzerland). Amplicons were purified using AMPure XP beads (BeckmanCoulter Inc., Fullerton, CA, USA). Indexing and sequencing (MiSeq platform, Illumina Inc., San Diego, CA, USA) were conducted by the Genomics Facility Basel (Basel, Switzerland). Library preparation and Illumina sequencing of 16S rRNA genes were performed by Microsynth AG (Balgach, Switzerland). 16S rRNA genes were amplified using the primers S-D-Bact-0341-b-S-17 (5′-CCTACGGGNGGCWGCAG-3′, 341F) and S-D-Bact-0785-a-A-21 (5′-GACTACHVGGGTATCTAATCC-3′, 805R; Herlemann et al., 2011). Anaerobic methanotrophic archaea (ANME) or related archaea were excluded from the analysis as previous research showed they are of minor importance in the Rotsee water column (Oswald et al., 2015).

16S rRNA gene and pmoA sequencing data were analyzed by the Genomic Diversity Centre (GDC, Zurich, Switzerland). Raw data were quality controlled using FastQC (v0.11.4) and MultiQC (v0.7). 16S rRNA gene low quality ends of reads were trimmed with PRINSEQ-lite (v0.20.4) and merged using usearch (v8.1.1812_i86linux64). pmoA reads were trimmed and merged using usearch (v9.2.64_i86linux 64) and FLASH (v1.2.11), respectively. Merged reads were primer-site trimmed by cutadapt v1.5 and v1.12, respectively. PRINSEQ-lite (v0.20.4) was used to filter and size-select the amplicons. Sequences were clustered into operational taxonomic units (OTUs) at a 97% similarity cut-off level using the uparse (16S rRNA gene) and an 86% similarity cut-off level using the unoise (pmoA) workflow with usearch (Wen et al., 2016). Taxonomic assignment of the representative sequences was set using the UTAX classifier together with the GreenGene database (May 2013) for 16S rRNA gene, and the SINTAX classifier with the pmoA database from (Wen et al., 2016). All sequence data have been deposited at the ENA database under the accession numbers PRJEB28460 (16S rRNA gene) and PRJEB28505 (pmoA).

As a basis for ecological interpretation and statistical analysis, we divided the water column into four zones with contrasting environmental conditions (Figure 1). The “oxic” zone ranged from the surface to depths where O₂ started to decrease and where CH₄ remained stable at low concentrations, but epilimnetic CH₄ sources might be present (Hofmann et al., 2010; Donis et al., 2017). This zone was followed by the “oxycline” zone, where O₂ was present but declined with depth, while CH₄ remained low. The lower boundary of the oxycline zone ended where CH₄ accumulation began. Below the oxycline, O₂ was below detection, but in Rotsee in-situ photosynthetic O₂ production still provides oxygen (Oswald et al., 2015). The “oxidation zone,” in which CH₄ was mainly consumed, was thus the zone from where CH₄ began to increase and either O₂ was present at measurable concentrations or light levels (PAR) allowed for photosynthesis. Therefore, this zone ended where PAR fell below detection limit. The dark bottom water was defined as the “anoxic” zone, where O₂ for aerobic CH₄ oxidation was lacking.

Statistical analysis was performed within the R statistical software environment (version 3.4.3; R Core Team, 2020) using the vegan, phyloseq, mixOmics, EcoSimR, igraph, and picante packages (Csárdi and Nepusz, 2006; Kembel et al., 2010; McMurdie and Holmes, 2013; Gotelli et al., 2015; Rohart et al., 2017; Oksanen et al., 2018). 16S rRNA gene (bacterial community) and pmoA (MOB) reads occurring less than 3 or 50 times in at least three samples, respectively, were removed from the OTU tables. To compare relative OTU abundances, reads were standardized to the mean sequencing depth. 16S rRNA gene data assigned to MOB species was excluded in statistical models incorporating the pmoA data set to avoid redundant correlations. Alpha diversity measures were calculated using default vegan functions (richness, Shannon). To obtain profiles for comparison and for co-linearity assessment, a number of missing values were imputed using the missMDA package [10 complete profiles for NO₂, NO₃, NH₄, SO₄, PO₄, DIC (all June 2013), Chl-a (September 2014 and 2015), δ¹³C-CH₄ (June 2013, September 2015), and eight missing values for pH in September 2014; star labeled plots in Supplementary Figure S1; Josse and Husson, 2016]. Except for pH, the imputed profiles were not incorporated in the statistical models but are discussed in the context of their co-linearity with non-imputed variables. Predicted values were tested by multiple‐ and over-imputation (Supplementary Figure S2).

A principal component analysis (PCA) with the environmental variables was performed to summarize the changes in physico-chemical variables along the depth gradients during temporal succession. Pearson correlations with Holmes-corrected p-values of physico-chemical parameters were calculated to assess the influence of single physico-chemical variables on specific MOB species. The mvabund package was used to test environmental, total bacterial and MOB community structure differences between different depth zones (Wang et al., 2012). To assess the environmental influence on the MOB community structuring, we applied canonical correlation analysis (CCA) based on selected environmental parameters (forward, backward, and both, incorporating variables that were chosen at least in two selection strategies). Variables with variance inflation factors (VIFs) > 10 were removed prior to model selection and significance was assessed by permutation tests (9,999 permutations).

Randomness in co-occurrence of MOB and total bacterial community was tested with C-score metric and quasiswap algorithm in a null model using the EcoSimR package (Gotelli and Ulrich, 2012; Gotelli et al., 2015). To determine if phylogenetically related species cluster within specific sites in different depth zones, we assessed the standardized (z-score) effect size (SES) of mean pairwise distances (MPD) and the mean nearest taxon distances (MNTD; Webb et al., 2002). SES were compared to a null model (“richness” null model, 999 randomizations of the phylogenetic tree) and differences between zones SES were assessed by ANOVA after checking for variance homogeneity and normality of the data distributions, followed by Tukey’s HSD. MPD is sensitive to differences in phylogenetically more distant taxa, whereas MNTD is sensitive to differences of phylogenetically more closely related taxa (i.e., tip of a phylogenetic tree). Non-metric multidimensional scaling (NMDS) of Bray-Curtis dissimilarity matrix of the square root transformed Wisconsin standardized OTU tables (bacteria and MOB) were determined. Overall similarity of MOB and the total bacterial community was assessed by Procrustes analysis of the first two dimensions of the respective NMDS plots. Significance of the Procrustes was estimated by a correlation-like statistic based on the squared m12 algorithm (9,999 permutations; Peres-Neto and Jackson, 2001).

To evaluate the inter-correlations of the total bacterial community, MOB and environmental variables explaining the depth zonation, we used the supervised N-integration discriminant analysis with DIABLO from the mixOmics package in R (Rohart et al., 2017). This analysis extracts complementary information from several data sets measured on the same samples but across different data type platforms (i.e., pmoA, 16S rRNA gene, and environmental variables) and aims to understand the interplay between the different levels of data that were measured. Optimal sparsity parameters were determined by computing M-fold cross-validation scores. The relative influence of inter-correlations of the selected environmental variables and bacterial community members on MOB at an association relevance level > 0.5 was assessed by partial Bray-Curtis dissimilarity-based redundancy analysis using Hellinger standardized 16S rRNA gene and pmoA data and scaled and centered environmental data. Variables having a VIF > 5 (environmental data) and VIF > 3 (16S rRNA gene) were stepwise removed and selected (forward, backward, and both) prior to analysis. The results of the DIABLO approach were analyzed by building relevance networks visualized with Cytoscape 3.7.2 using the EClerize layout. Non-randomness of the networks was tested by comparing the network to 10,000 random Erdös-Réyni networks with similar numbers of edges and nodes (Ju et al., 2014; Weiss et al., 2016). The network GML files are available as Supplementary Materials. See also section Supplementary Material for additional statistical analysis and plots.

The defined depth zones over all four field campaigns were consistently different from each other in terms of their overall physico-chemical properties (mvabund: Likelihood Ratio Test (LRT) = 778.8, p < 0.001; Figure 1), and sampling time contributed considerably to the variability in the dataset (see ordination of samples based on their physicochemical characteristics in Figure 2). The annual T-driven stratification during summer months resulted in a narrowing of the oxycline and the oxidation zone as the season progressed (see T(°C) profiles in Supplementary Figure S1). The surface water was always well-oxygenated and O₂ concentrations fell below detection limit usually in the upper half of the oxidation zone (Figure 1). The anoxic zone was substantially enriched with CH₄, while only small amounts of CH₄ (0.12–1.07 μM) were detected in the oxic zone, which, however, was still oversaturated relative to the atmosphere (Schubert et al., 2010). In September 2015, the oxycline was combined with the oxidation zone as CH₄ profiles already slightly increased, where O₂ decreased. In June 2013, light was detected almost to the sediment, therefore an anoxic zone was not defined. Vertical profiles of additionally measured physico-chemical parameters are summarized in Supplementary Figure S1. Turbidity (Turb) maxima were mostly congruent with Chl-a throughout the water column as profiled in June and August 2013. The anoxic zone exhibited substantial concentrations of TDN, PO₄, DIC, and reduced substances (Fe_DGT/Diss, Mn_DGT/Diss, S_Tot, and NH₄). Cu_DGT, Cu_Diss, and Cu_Tot were found to be highest at the lake’s surface and decreased strongly in the lower oxic zone and in the oxycline, whereas Cu_Part concentrations usually peaked within the oxidation zone. Some variables exhibited pronounced co-linearity (see vectors in PCA Figure 2; Supplementary Figure S2).

An average of 41,951 reads per sample were assigned to 1,829 unique bacterial OTUs after filtering. Alpha diversity over all campaigns increased from around 400 OTUs in the surface water to approximately 1,000 OTUs close to the sediment (see alpha diversity plots in Supplementary Figure S3). Sixteen OTUs were assigned to proteobacterial MOB and five OTUs to potential verrucomicrobial MOB. Together they represented 1.07% of all filtered reads. These OTUs were removed in analysis for co-occurrence that incorporated both 16S rRNA gene and pmoA data to avoid trivial correlations. Actinobacteria, Bacteroidetes, non-MOB Proteobacteria, and Verrucomicrobia dominated the communities during all campaigns (Figure 3). Actinobacteria and Proteobacteria abundance tended to decrease or increase with depth, respectively. Cyanobacteria were abundant in both the oxic zone and oxycline in September 2015 and were present in lower proportions in the other three campaigns. Firmicutes were detected in the oxidation and anoxic zone. OD1 were highly abundant in the oxidation and anoxic zone in August 2013 and in September 2014 and 2015. Planctomycetes inhabited the oxic zone and oxycline mainly in August 2013. Microbial community structures were different between the sampling dates (mvabund: LRT = 29,741, p < 0.001) and not randomly distributed among depths over all campaigns (C-score metric for randomness in co-occurrence = 14.16, p < 0.001, SES = 24.75). Bacterial communities were further structured along the depth gradient during all campaigns and were significantly different between the depth zones (mvabund: LRT = 45,650, p < 0.001; see NMDS ordination in Supplementary Figure S4).

We obtained an average of 89,816 pmoA reads per sample, which resulted in 121 OTUs after removing sparse OTUs (i.e., less than 50 reads in at least three samples). Alpha diversity increased with depth from 21 OTUs in the surface waters to 110 OTUs in deeper waters (Supplementary Figure S3). Sequences from type I (Gammaproteobacteria) and type II (Alphaproteobacteria) MOB were identified (Figure 3). Type I MOB were assigned to Methylobacter, Methylomonas, Methylosoma, and various environmental clusters (type Ia and type Ib), whereas type II MOB comprised only Methylocystis. MOB communities were structured along the depth gradient and differed between the depth zones and campaigns (mvabund: LRT = 3,134/2,900, p < 0.001; see NMDS ordination in Supplementary Figure S5). June 2013 was dominated by Methylobacter and type Ia (herein, we refer to type Ia excluding Methylobacter, Methylomonas, and Methylosoma) predominantly inhabiting the lower oxycline and the oxidation zone (Figure 3). Methylobacter was abundant in August 2013 in the oxidation and anoxic zone, whereas type Ia occurred from the oxycline on downward. In September 2014 and 2015, Methylobacter was also detected within anoxic waters. Methylocystis was found mainly in the oxic part of the lake, with the highest relative abundance in September 2015. Methylomonas and Methylosoma were most abundant in the CH₄ oxidation zone of September 2014. Type Ib MOB were detected in low numbers in August 2013 in the lower part of the oxic zone and could also be found more dispersed within the oxic zone and the oxycline in September 2014 and 2015.

Statistical testing confirmed that MOB communities were not randomly distributed within the water column (C-score = 118.93, p < 0.001, SES = 98.4) and showed phylogenetic relatedness higher than expected by chance in each sample within the different zones (SES of MPD < 0, SES of MNTD < 0, SES are shown in Supplementary Figure S6). Phylogenetic relatedness at higher node levels changed between specific samples in the oxycline compared to the anoxic zone [ANOVA F(3,66) = 3.70, p < 0.05, Tukey’s HSD < 0.01], whereas phylogenetic relatedness at the tree-tip level (i.e., lower node levels) was highest within the oxidation and anoxic zone [ANOVA F(3,66) = 7.17, p < 0.001, Tukey’s HSD < 0.01]. These statistics indicate that MOB clades that are phylogenetically highly similar coexist within specific sites in the oxycline, whereas phylogenetically less similar MOB tend to be mutually exclusive. Inversely, sites within the anoxic zone have a narrower phylogenetic structuring but highly similar MOB tend to be mutually exclusive. Broad phylogenetic correlation with environmental structuring is evident, but some MOB show environmental preferences that are distinct from their close phylogenetic relatives (see color coded bar for genus in Figure 4).

A parsimonious CCA based on selected environmental variables explained a total of 69% of the MOB community structure over all campaigns with 42% of variation being explained on the first two axes (Figure 5). CH₄ and O₂ were the expected strong antipodal drivers. Light, Cu_Diss, pH, Fe_DGT, and TDN also contributed to the first canonical axis. Procrustes analysis showed statistically significant similarity between structuring of MOB and the rest of the bacterial community along the depth gradients over all sampling dates (Procrustes correlation = 0.51, p < 0.001). This could indicate either biological interactions between MOB and other microbes or similar niche preferences.

Network analysis provides a framework to study associations between several classes of variables (Barberán et al., 2012; Comte et al., 2016). We constructed a relevance network to analyze the connection of MOB and the bacterial community composition (pmoA data and 16S rRNA gene data excluding 16S rRNA gene MOB OTUs from the latter) and physico-chemical variables (Figure 6A, for network characteristics see Supplementary Table S1). Fifty-nine MOB, 271 bacterial OTUs, and 17 physico-chemical variables formed distinct networks, with sub-networks that conformed to the depth zonation of Rotsee. This confirms that the originally hypothesized zonation can be broadly reconstructed from the dataset. The network architecture showed positive inter-correlations of Methylocystis with O₂, Cu_Diss, light, pH, T, and Zn_Diss in the surface waters (oxic zone). Type Ia MOB in the oxycline were connected to Mn_Part and Cu_DGT, and with Cu_Part and Turb in the oxidation zone. Within greater depths (anoxic zone), Methylobacter was linked to a larger set of bacteria and physico-chemical variables (Mn_Diss, CH₄, Fe_Diss, TDN, S_Tot, Fe_DGT, and Mn_DGT). Connections of MOB with bacterial OTU were common in all depth zone subclusters and are summarized in Figure 6A (i.e., number of connections of specific phylum within the network).

We used partial redundancy analysis to distinguish the relative contributions of physico-chemical variables and bacterial interactions on the variance of MOB abundance. According to this analysis, out of a total of 84% of explained variance in MOB occurrence, 22% can be explained exclusively by bacterial interactions, but only 2% exclusively by environmental drivers (Figure 6B). Interestingly, most of the explained variance is shared by physico-chemical variables and the bacterial community (60%). Specific bacterial OTUs affiliated with MOB in the network analysis are illustrated in the Cytoscape file (Supplementary File “Relevance network.cys”).

The availability of Cu as a co-factor of pMMO’s active site can restrict enzymatic activity in MOB communities and thus limit their growth (Semrau et al., 2010, 2018). Previous work on the role of Cu for MOB in Rotsee has already established a likely role of Cu scarcity and competition for Cu in the lake (Guggenheim et al., 2019). The present work provides additional evidence for Cu as an important factor controlling MOB community assembly: we observed positive correlations of MOB with changing importance of Cu-species along the depth gradient of Rotsee. Cu_Diss correlated with Methylomonas and type Ib OTUs, but most strongly with Methylocystis in the surface water (Figures 4, 5). These genera are able to use Cu acquisition mechanisms based on complexing agents to deal with low bioavailable Cu supply conditions (Ul-Haque et al., 2015). Such auxiliary peptides could mobilize Cu from the non-bioavailable Cu_Diss fraction (Dassama et al., 2017; Kenney and Rosenzweig, 2018). Additionally, most members of this taxon possess a high CH₄ affinity copper-dependent pMMO isozyme, which would support oxidizing CH₄ at the sub-micromolar CH₄ concentrations prevalent in the surface water (Baani and Liesack, 2008; Reis et al., 2020). Methylocystis is often found in warmer waters (Borrel et al., 2011; Tsutsumi et al., 2011) and it is suggested that high T selects for type II over type I MOB (Sundh et al., 2005). Indeed, the warmer surface water of Rotsee seems to favor the presence of Methylocystis. It is thus likely that Methylocystis contribute to CH₄ consumption in the oxic zone.

During late stratified periods, when CH₄ has accumulated in the hypolimnion of Rotsee, the highest MOB abundance was found at the lower end of the oxycline and in the oxidation zone (Figure 3). Under these conditions Cu_DGT and Cu_Part correlated mostly with type Ia and Methylobacter OTUs (Figures 4, 5). It is thought that Cu_DGT is the bioavailable Cu fraction, whereas Cu_Part mainly represents Cu incorporated into biomass (Guggenheim et al., 2019). Cu_DGT concentrations decrease substantially from the oxycline toward the oxidation zone (Supplementary Figure S1). In situ studies focusing on the influence of Cu on MOB community structures are rare, but it has been elucidated that MOB able to express pMMO thrive even under very low levels of bioavailable Cu (< 50 nM; Cantera et al., 2016). In addition, as mentioned before, certain MOB species possess special Cu uptake mechanisms to increase the bioavailable Cu fraction.

Previous studies also reported the occurrence of Methylobacter species in the anoxic zones of stratified lakes (Biderre-Petit et al., 2011; Blees et al., 2014; Milucka et al., 2015; van Grinsven et al., 2020). Bioavailable and dissolved Fe (Fe_DGT and Fe_Diss) correlated with most Methylobacter species (Figures 4, 5), although cultivated representatives are not able to express the iron-dependent sMMO enzyme (Knief, 2015) and sMMO appears to be rare in Rotsee (Guggenheim et al., 2019). However, it is possible that these MOB rely on Fe for other enzymatic pathways (e.g., formate dehydrogenase; Glass and Orphan, 2012). Furthermore, it has been suggested that the mechanism of Fe-coupled anaerobic CH₄ oxidation is accomplished by a complex microbe-mineral reaction network in which both, MOB and iron-reducing organisms (bacteria and archaea) are directly and indirectly involved (Bar-Or et al., 2017; Cabrol et al., 2020). For example, besides methanogens being able to produce CH₄, they are additionally involved in reducing Fe-oxides at high CH₄ concentrations leading to intermediates, which are required by MOB for CH₄ oxidation (Bar-Or et al., 2017). Mn_DGT and Mn_Diss also correlated with Methylobacter OTUs within the anoxic zone (Figures 4, 5). It has been proposed that an alternative anaerobic CH₄ oxidation lifestyle proceeding via Fe(III) or Mn(IV) reduction could be relevant in CH₄-rich anoxic zones of lakes, which could explain the increase in Fe_Diss and Mn_Diss in the deeper waters of Rotsee and its correlation with Methylobacter (Crowe et al., 2008; Oswald et al., 2016a; Supplementary Figure S1).

Interestingly, Zn_Diss was also suggested as a significant driver by our analysis and was not co-linear with other environmental variables (Figures 4, 5; Supplementary Figure S2). There was a positive correlation of Methylocystis with Zn_Diss although Zn potentially inhibits pMMO activity (Sirajuddin et al., 2014). However, the highest measured Zn_Diss concentrations in Rotsee were ~100 times lower than those tested experimentally by Sirajuddin et al. The close network connectivity within the oxic zone between Methylocystis and other bacteria (i.e., linkage to Actinobacteria, Bacteria, Proteobacteria, and Bacteroidetes; Figure 6B) might indicate indirect effects mediated by Zn. For example, the bacterial OTU1 could be assigned to the order Candidatus Nanopelagicus within Actinobacteria isolated from Lake Zurich (Neuenschwander et al., 2018; see also the cytoscape file). The isolate showed a reduced genome and thus might have strong metabolic dependencies on co-occurring bacteria (i.e., MOB) for lost metabolic functions that have to be provided by functional leakage (Morris et al., 2012). As Actinobacteria exhibit Zn concentration linked gene regulation mechanisms, this might mirrors a reverse effect of Methylocystis on this Candidatus Nanopelagicus (Choi et al., 2017), mediated by, i.e., the release of riboflavin, nicotinamide, and thiamine of the Candidatus Nanopelagicus. Anyway, further work will be necessary to confirm the role of Zn and to investigate potential mechanisms linking its availability to MOB dynamics.

The highly significant contribution of total dissolved nitrogen (TDN) as an explanatory variable for MOB abundance indicates links between MOB and nitrogen availability (Figures 4, 5). In particular, TDN contributed to the position of Methylobacter in the relevance network (Figure 6A). The observed gradient of TDN strongly correlated with NH₄ and exhibited similar concentration ranges suggesting that NH₄ contributes the largest part of TDN (Supplementary Figures S1, S2). NH₄ is a central nutrient in aquatic systems, hence a positive correlation between NH₄ and MOB could be due to the fact that MOB assimilate NH₄ for growth. Previously, a laboratory study with littoral wetland samples from a boreal lake has demonstrated that nitrogen load in form of NH₄NO₃ changed the MOB community structure and favored activity of type I MOB, particularly Methylobacter cells (Siljanen et al., 2012). The evolutionary linkage of the genetic sequence of pmoA and amoA, which encodes for the ammonia monooxygenase (AMO), endows most MOB, especially type I MOB, with the ability to oxidize NH₄ through the pMMO enzyme (Knief, 2015; Khadka et al., 2018), which may be another explanation for the link between MOB and NH₄.

CH₄ oxidation can also be coupled directly to the nitrogen cycle. Anaerobic oxidation of CH₄ via NO₂-/NO₃-reduction (n-damo) has normally been attributed to Methylomirabilis species of the NC10 phylum, which according to previous data do not appear in Rotsee (Mayr et al., 2019). However, recent whole genome and environmental metagenome analysis have revealed that various assimilatory and dissimilatory nitrogen reduction genes, such as those encoding NO₂‐ and NO₃-reductases, are also found among gammaproteobacterial MOB, especially Methylobacter and Methylomonas species (Chistoserdova, 2015; Zhu et al., 2016). This would suggest that NO₂ and NO₃ can be used by MOB to oxidize CH₄. However, NO₂ and NO₃ were not detected below the oxidation zone in Rotsee (Supplementary Figure S1).

Spatio-temporal shifts in the MOB community were highly congruent with the changing composition of the total bacterial community (OTU level, Procrustes 51%). The restrictive and congruent localization of MOB and bacterial OTUs throughout the stratified periods suggests that biological interactions between MOB and other bacteria need to be analyzed in more detail.

In our network analysis we focused exclusively on positive MOB co-occurrence patterns with environmental variables and the bacterial community. This means that observed associations can arise from a number of causes: species co-occurrence can be driven by direct biological interactions such as synergism (i.e., cooperation) but also by similar niche preferences (neutral effects; Faust and Raes, 2012; Ho et al., 2016). Additionally, positive feedback loops between MOB and the total bacterial community interacting with their environment may also result in co-occurrence. It is noteworthy that positive interactions with bacteria alone (excluding impact of the environment) explained a relevant part of the MOB occurrence (22%; Figure 6B). This indicates that at least some of the correlations are through actual microbial interactions and do not just reflect overlapping niche preferences. Anyway, it is possible that these interactions are mediated through an effect chain of non-measured variables, i.e., by the influence of trophic levels ultimately influencing MOB dynamics.

However, the majority of variance in MOB abundance within the network was explained jointly by physico-chemical parameters and the bacterial community (60%; Figure 6B), while according to CCA environmental parameters explained 69% of MOB abundance. While this underscores the importance of physico-chemical parameters as drivers of the bacterial and MOB community structure, positive correlations of this type could still involve biological interactions. Physico-chemical variables may, e.g., affect the bacterial community structure, which subsequently influences the MOB community by synergistic effects or vice versa. The available data do not allow us to distinguish between these possibilities. It is further possible that unmeasured environmental variables contribute to shaping of the MOB community structure. However, it is noteworthy that bacterial and physico-chemical inter-connectivities showed high complexity in the surface layer and the hypolimnion of the lake (i.e., number of connections in Figure 6A), while inter-connectivities within the oxycline and oxidation zones were less distinct.

Hereafter, we suggest some examples of possible positive bacterial interaction mechanisms taking place in the water column of Rotsee. Within the oxidation zone, part of the type Ia and Methylobacter communities correlated with Turbidity (Turb), a possible indicator for primary producers as Turb showed some overlapping patterns with Chl-a depth profiles (Supplementary Figure S2). MOB can form mutualistic interactions with oxygenic phototrophs in light penetrated anoxic layers enabling CH₄ oxidation, while potentially providing carbon dioxide in return (Milucka et al., 2015). It seems likely that these MOB clusters are the main CH₄ consumers in Rotsee as the main aerobic CH₄ oxidation process in Rotsee during stratification happens within the oxycline and oxidation zone and might be predominantly coupled to oxygenic primary production (Oswald et al., 2015; Brand et al., 2016). There is also evidence that CH₄ oxidation by MOB under O₂ limitation is indirectly connected to the reduction of alternative terminal electron acceptors by other organisms in the anoxic waters of the lake. NO₂ and NO₃ were not measured in the anoxic waters of Rotsee, however, NO₂ could have been produced by Fe-/Mn-dependent anaerobic NH₄-oxidizing bacteria, and readily used by gammaproteobacterial MOB for oxidizing CH₄ (Ferousi et al., 2017; Kuypers et al., 2018). Due to the previously mentioned homology between pMMO and AMO, MOB and nitrifying bacteria are capable of both, CH₄ and NH₄ oxidation. Hence, the coexistence of MOB and nitrifying bacteria under anoxic conditions could be explained on the basis of similar niche preferences (Costa et al., 2019). Co-occurrence of MOB in the same areas as non-methanotrophic methylotrophs might be beneficial for MOB. Since high methanol concentrations inhibit MOB performance, co-occurring methylotrophs, which are able to assimilate methanol, remove this compound to levels that enable MOB to thrive (Denfeld et al., 2016). A recent in-situ study shows that non-methanotrophic methylotrophs induce a change in the expression of MOB MDHs via putative secretory compounds leading to an increased loss of methanol, which is readily taken up by the methylotrophs (He et al., 2012b; Biderre-Petit et al., 2019). For example, we found highly connected methylotrophs belonging to the family Methylophilaceae and the Verrucomicrobia subdivision 6 in the oxic zone in the relevance network. In the anoxic zone, we found members of the S-BQ2-57 soil group belonging to Verrucomicrobia that showed high connectivity within the network. Methanol also plays a significant role between the coupling of aerobic CH₄ oxidation and denitrification by the cooperation between MOB and denitrifying bacteria. Organic metabolites (i.e., methanol, citrate, acetate, formaldehyde, and formate) released by MOB could serve as electron donors for denitrification, where methanol is thermodynamically considered as the ideal trophic link (Kalyuhznaya et al., 2009; Zhu et al., 2016).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.