A metagenomic insight into freshwater methane-utilizing communities and evidence for cooperation between the Methylococcaceae and the Methylophilaceae

Experimental setup, sample collection and stable isotope probing

The schematic of the experiments conducted is depicted in Fig. 1. Sediment samples were collected on March 3, 2009, from a 63 m deep station in Lake Washington, Seattle, Washington (47.038075’ N, 122.015993’ W) using a box core that allowed collection of undisturbed sediment. Samples were transported to the laboratory on ice and immediately used to set up microcosms. In order to assess populations active in methane oxidation under different oxygen tensions and to test for their potential dependence on the presence of nitrate we set up microcosm incubations as follows. One microcosm was incubated in an atmosphere of 50% ¹³C-labeled methane (99 atom % ¹³C, Sigma-Aldrich) and 50% ambient air, to assess the populations active in aerobic methane oxidation (the +O₂-NO₃^- condition); the second microcosm was incubated in an atmosphere of 50% ¹³C-methane and 50% ambient air, in the presence of 10 mM KNO₃, to assess the populations active in aerobic methane oxidation positively responding to the presence of nitrate (the +O₂+NO₃^- condition); the third microcosm was incubated in an atmosphere of 50% ¹³C-methane and 50% N₂, to assess the populations active in microaerobic methane oxidation (the -O₂-NO₃^- condition); the fourth microcosm was incubated in an atmosphere of 50% ¹³C-methane and 50% N₂, in the presence of 10 mM KNO₃, to assess the populations active in microaerobic methane oxidation positively responding to the presence of nitrate (the -O₂+NO₃^- condition). Each microcosm contained 50 ml of the top layer (1 cm) of the sediment and 50 ml of Lake Washington water filtered through 0.22 µm filters (Millipore). Samples were placed into 250 ml glass vials (6 vials per experiment, the contents of which were mixed before DNA extraction) and these were sealed with rubber stoppers. The duration of the incubation time for each microcosm was determined empirically by observing the formation of a heavy, ¹³C-enriched DNA fraction. It took 10 days for heavy DNA band to appear in the +O₂+NO₃^- microcosm, compared to 15 days for the +O₂-NO₃^- microcosm, suggesting that the methane-consuming community was stimulated by nitrate. It took much longer for heavy DNA band to appear in the microaerobic microcosms (20 and 30 days, respectively, for nitrate-amended and nitrate-free conditions; Fig. 1). These data suggest that nitrate also had a positive effect on methane consumption by the microbial community in microaerobic conditions.

Community DNA was extracted as described previously (Beck et al., 2011) with one modification as follows: DNA samples were subjected to an additional round of purification using UltraClean^® Mega Soil DNA Isolation Kit (MOBIO). The heavy (¹³C-enriched) fractions of DNA were separated from the light (¹²C) fractions by CsCl-ethidium bromide density gradient ultracentrifugation, visualized under UV (Fig. 1), and collected and purified following standard procedures as previously described (Neufeld et al., 2007). DNA from an untreated sediment sample was extracted using the same protocol. This sample was not subjected to density gradient centrifugation.

DNA sequencing and assembly

Five shotgun libraries were constructed, one from each microcosm. The libraries were sequenced using the 454 sequencing technology at the Joint Genome Institute (JGI) Production Genomics Facility. A total of 5,241,266 reads comprising 1.67 gigabases (Gb) of sequence were generated. These were assembled using the Newbler assembler. Assembly statistics are shown in Table 1.

Pyrotag sequencing

We used 454 pyrotag sequencing of the V8 hypervariable region of the 16S rRNA gene to determine the compositions of microbial communities in the four ¹³C-enriched metagenomes using a computational pipeline PyroTagger, as previously described (Kunin & Hugenholtz, 2010), and we compared these to the community composition in the original (unamended) lake sediment sample. The number of total pyrotag sequences generated per microcosm is shown in Table 1. The PyroTagger assignments were used to calculate the Shannon-Weaver diversity index and the expected number of unique species rarefied to a constant value across all samples (Table 1).

Metagenome analysis

The draft quality assemblies were processed using the IMG/M pipeline (Markowitz et al., 2012) and the outputs of these automated analyses were manually verified and used for metagenome profiling.

Single gene taxonomy

To classify the 16S rRNA gene sequences in the metagenomes, the 16S genes identified by JGI pipeline were aligned against the Ribosomal Database Project (RDP) version 10, release 26. The top scoring alignment for each sequence was used to assign taxonomy, based on the annotations in the RDP. Genes were classified to the family level.

To detect pmoB and mmoX genes, respective peptides representing Methylococcaceae and Methylocystaceae were obtained from public databases and used as queries in pBLAST analyses against each dataset, and all matches were recorded. These queries pick each other and also the NC10 sequences. All matches were compared with sequences in public databases using pBLAST. Genus level assignments were made using an 80% identity cutoff level, and family level assignments were made at 70% identity cutoff, based on prior knowledge on gene diversity within each family (Chistoserdova & Lidstrom, 2013).

To detect fae genes, respective peptides representing Methylococcaceae, Methylocystaceae, and Methylophilaceae were obtained from public databases and used as queries in pBLAST analyses against each dataset, and all matches were recorded. These queries pick each other, other proteobacterial sequences as well as NC10, planctomycete and unclassified sequences. Genus and family level assignments were done as above.

For the nitrogen metabolism genes, sequences were selected by their annotation and aligned against the NCBI’s nonredundant database (nt). The top scoring hit for each sequence was saved. The NCBI taxon ID was extracted and the NCBI taxonomy database queried to collect the “scientific name” of each hit. The count of hits for each “scientific name” were reported.  The following annotation search terms were used: “itrate reductase” for NO₃-reductase, “itrite reductase” for NO₂-reductase, “itric[-]oxide reductase” for NO reductase, “itrous[-]oxide reductase for N₂O reductase, and “itrogenase” for nitrogenase. The first letter was omitted to avoid conflicts with upper and lower case letters and in the case of NO and N₂O reductases, both annotations (with and without “-”) were accepted.

To classify mxaF and xoxF genes belonging to specific methylotroph families, high stringency pBLAST analyses were applied using 90% identity level cutoff at the protein level and employing publically available sequences of the respective enzymes representing the respective taxonomic groups.

Pyrotag profiling of community DNA shows enrichment for Methylococcaceae and Methylophilaceae sequences in aerobic microcosms

As expected from prior analyses (Kalyuzhnaya et al., 2008), the community in the unamended sample revealed high complexity, being represented by a total of 1,486 sequence clusters (97% sequence identity; Kunin & Hugenholtz, 2010). The community was dominated by Proteobacteria (33.3%), of which phylotypes of the Methylococcaceae family that represents one class (called type I) of methane oxidizing bacteria were most prominently present (10% of all sequences). The second most dominant group was represented by chloroplast sequences (21.9%). Other prominently present phyla were Bacteroidetes (10.5%), Acidobacteria (7.2%) and Chloroflexi (4.0%; Fig. 2A). The phylogenetic complexity of the ¹³C-enriched metagenomes representing three of the enrichment conditions (the +O₂-NO₃^- condition, the +O₂+NO₃^- condition, and the -O₂-NO₃^- condition) was significantly reduced compared to the non-enriched community (313, 709 and 561 sequence clusters, respectively), and different shifts in phyla distribution occurred in these communities. The proportion of proteobacterial sequences increased in the aerobic communities (to 49.7% and to 82.0%, respectively, in +O₂-NO₃^- and in +O₂+NO₃^- conditions; Fig. 2A). In both cases phylotypes classified as the Methylococcaceae and the Methylophilaceae were most prominently present (10.9% and 14.9% of all sequences in the +O₂-NO₃^- condition and 56.5% and 10.5% of all sequences in the +O₂+NO₃^- condition, respectively; Fig. 2B). The proportion of Proteobacteria decreased (to 13.2%) in the -O₂-NO₃^- condition, the dominant phylotypes being the chloroplast and the Bacteroidetes phylotypes (34.0% and 9.5% of total sequences, respectively), while the Methylococcaceae and the Methylophilaceae phylotypes constituted only a minor fraction of all sequences (0.03% and 0.9%, respectively). The proportion of proteobacterial phylotypes and the overall make up of the community of the -O₂+NO₃^- condition, at the phylum level, resembled that of the non-enriched community (Fig. 2B). However, the Methylococcaceae phylotypes were less represented (4.0% of total sequences) while other proteobacterial phylotypes (such as Burkholderiales) were more represented than in the non-enriched community. Phylotypes representing other bona fide methylotroph taxa, including members of the Methylocystaceae and Bejerinckiaceae that represent the second class (called type II) of methane oxidizers and members of the recently described NC10 phylum implicated in anaerobic methane oxidation linked to denitrification (Wu et al., 2011) were detected in all five pyrotag libraries. However, their proportions were very small, not exceeding 0.5% of the total community in each case.

Genome recruitment further highlights the dominant presence of Methylococcaceae and Methylophilaceae species in aerobic microcosms

A total of 1,362,213,455 base pairs (1.36 Gb) of assembled sequence were generated. Genes were called using the standard JGI IMG/M pipeline (Markowitz et al., 2012), and each gene was taxonomically classified by using its best BLAST hit in the current genomic database employed by the IMG/M interface (3811 bacterial genomes, 163 archaeal genomes, 177 eukaryotic genomes and 2803 viral genomes; sequencing and assembly statistics are shown in Table 1). For taxonomic assignments, we considered separately all protein coding genes and genes classified at the 60% and the 90% cutoff levels (protein level classification). While few genes were classified at the 90% cutoff level (2.44 to 17.56% of the total, dependent on the microcosm) these provide very robust proxies for the organisms represented in the metagenomes, especially given the fact that genomes of key model methylotrophs were parts of the database used for comparisons, including genomes originating from Lake Washington (Lapidus et al., 2011; Kittichotirat et al., 2011 and unpublished). Thus we mostly relied on the 90% cutoff classification for the confident sequence assignments, while realizing that these provide the lowest estimate for the presence of a specific phylum. At the 60% cutoff level, between 28.9% and 42.5% of the total genes could be classified, allowing for encompassing species not represented in the databases by very closely related models. In either case, most of the genes in each metagenome were taxonomically classified as proteobacterial, with the absolute majority matching to beta- and gammaproteobacteria (up to 45.4% and up to 50.1% of genes at 90% cutoff level, respectively). The next most abundant group were alphaproteobacteria (up to 12.4% of genes at 90% cutoff level), excepting the -O₂-NO₂^- metagenome, in which, along with gamma- and betaproteobacteria, Bacteroidetes and deltaproteobacteria sequences were prominent. Remarkably, most of the gammaproteobacterial sequences (up to 84.2/98.4% at the 60/90% cutoff levels, dependent on the microcosm) were classified as belonging to the Methylococcaceae family. Among betaproteobacterial sequences, the largest proportion of sequences in the aerobic microcosms (up to 53.6/88.4%) were classified as belonging to the Methylophilaceae family, while in the microaerobic microcosms betaproteobacterial sequences were distributed among a number of dominant families, which included Methylophilaceae, Comamonadaceae, Rhodocyclaceae and Burkholderiaceae. Among the alphaproteobacterial sequences, those classified as belonging to the Methylocystaceae family only constituted the highest proportion (48.6/81.3%) in the +O₂-NO₃^- microcosm, but they were also abundant in the +O₂+NO₃^- microcosm (13.0/30.1% of total alphaproteobacterial sequences). In the remaining datasets, few alphaproteobacterial sequences were classified as Methylocystaceae, with Bradyrhizobiaceae sequences being most abundant. Similarly to the results of pyrotag analysis, sequences of other bona fide methylotrophs (such as Methylobacteriaceae, Bejerinckiaceae, Xanthobacteriaceae, Hyphomicrobiaceae) were identified in each metagenome. However, the proportions of the genes assigned to each of these families were low, suggesting that these species were minor members of the functional communities investigated.

Sequences assigned to each of the major methylotroph families, Methylococcaceae, Methylophilaceae and Methylocystaceae, were further classified at the genus level, by matching each gene classified to the family level to the available genomic scaffolds representing each family, at 60 and 90% cutoff values. In these analyses, the Methylococcaceae family was represented by a total of six genomes, of Methylobacter tundripaludum (Svenning et al., 2011), Methylomonas sp. (unpublished), Methylomicrobium album (unpublished), Methylocaldum szegediense (unpublished), Methylococcus capsulatus (Ward et al., 2004) and Crenothrix polyspora (unpublished). In the case of the latter, we find the claimed affiliation of this strain within a proposed new family of Crenothrichaceae (Stoecker et al., 2006) invalid, based on high sequence similarity with the members of Methylococcaceae (at least 95% at 16S rRNA gene level; (Iguchi, Yurimoto & Sakai, 2011), and thus we consider this organism as part of this family. However, the matters are complicated further by the fact that the sequence of the yet uncultivated C. polyspora originated from a highly enriched but not axenic culture that may contain other representatives of Methylococcaceae. Of the six genomes, only one represented an organism originating from Lake Washington, the Methylomonas sp. strain (unpublished data). The Methylocystaceae family was represented by only two genomes, of Methylosinus trichosporium (Stein et al., 2010) and Methylocystis sp. (Stein & Klotz, 2011), none originating from Lake Washington. The Methylophilaceae family was represented by eight genomes, of Methylotenera mobilis, Methylotenera versatilis, Methylovorus glucosotrophus (Lapidus et al., 2011), Methylophilus sp., unclassified Methylophilaceae (unpublished), Methylobacillus flagellatus (Chistoserdova et al., 2007), and unclassified Methylophilaceae strains HTCC8121 (Giovannoni et al., 2008) and KB13 (unpublished). Of these, the first five strains originated from Lake Washington (Kalyuzhnaya et al., 2006; Kalyuzhnaya et al., 2012 and unpublished) and the latter two were marine Methylophilaceae with unusually small genomes (Giovannoni et al., 2008). All the genomes mentioned here as “unpublished” have been sequenced by the JGI and are publically available through the IMG/M interface (http://img.jgi.doe.gov). From matching to this limited number of genomic scaffolds, the Methylobacter species appeared to be the dominant type among the Methylococcaceae representatives, the Methylocystis types appeared to dominate over the Methylosinus types, and Methylotenera species appeared to be dominant among the Methylophilaceae. However, specific dynamics could be observed at the genus level in response to different stimuli (Table 2). The proportion of the Methylomonas and Methylomicrobium types increased in response to the addition of methane, especially pronounced in the +O₂+NO₃^- condition, while the proportion of the “Crenothrix” types decreased compared to unamended sediment. The Methylobacter types were most dominant in the -O₂+NO₃^- condition. Dynamics were also clearly seen among the Methylocystaceae, with the proportion of Methylosinus sequences upshifting in response to methane in aerobic conditions and downshifting in response to nitrate, with respect to Methylocystis sequences. The relative proportion of Methylotenera sequences among the Methylophilaceae populations upshifted significantly in response to nitrate in aerobic conditions. We were also able to distinguish between two different species within Methylotenera, M. mobilis vs. M. versatilis, based on their significant divergence at the genomic level (Lapidus et al., 2011). While the proportion of M. versatilis-like sequences increased in the +O₂+NO₃^-condition, their proportion decreased in the -O₂+NO₃^- condition, being replaced by the M. mobilis-like sequences. Dynamics could also be seen among sequences classified as other Methylophilaceae. Unsurprisingly, virtually no sequences were matched to the scaffolds representing marine Methylophilaceae.

Taxonomic profiling demonstrated a good correlation between the populations of the Methylococcaceae and the Methylophilaceae in both aerobic conditions (Figs. 2 and 3). Both were more abundant in the +O₂+NO₃^- microcosm and somewhat less abundant in the +O₂-NO₃^- microcosm, suggesting that both preferred higher nitrate concentrations for C1 metabolism. The low oxygen tension conditions selected against all methylotroph species. However, the Methylococcaceae still represented the majority of gammaproteobacterial sequences at the 90% cutoff level (Fig. 2D). Overall, taxonomic profiling of metagenomes correlated well with pyrotag-based profiling both suggesting that Methylococcaceae and Methylophilaceae efficiently consumed the ¹³C label from methane in aerobic conditions while the label distributed more evenly among multiple phyla in microaerobic conditions. Results of ordination analysis of dissimilarity of the five communities are shown in Supplemental Figure 1.

Single gene-based taxonomic profiling supports data from whole-metagenome profiling

16S rRNA gene profiling in each microcosm was carried out (Tables 1 and 3; Supplemental Table 1). For the metagenomes with significant sequence sampling (Table 1), the distribution of 16S rRNA genes among major phyla matched well those determined by the pyrotag sequencing approach, with some differences such as the reduced proportion of chloroplast sequences (data not shown), which is likely due to the low diversity of chloroplast sequences. Analysis of 16S rRNA gene sequences revealed that only in the aerobic microcosms did methylotroph sequences make up a significant fraction of total 16S rRNA gene sequences (26.3 to 31.8%, respectively, in the +O₂-NO₃^- and the +O₂+NO₃^- conditions; Table 1). In the microaerobic conditions and in the unamended sample, the methylotroph 16S rRNA gene fraction made up approximately 4% of the total 16S rRNA sequences. The methylotroph sequences represented three major families, Methylococcaceae, Methylocystaceae and Methylophilaceae. Within each family, a variety of phylotypes were detected suggesting complex community composition within each class. While the Methylococcaceae sequences were most numerous in each microcosm (50 to 100% of the methylotroph 16S rRNA sequences), including the unamended sample, shifts in phylotype composition occurred in response to different incubation conditions. The +O₂+NO₃^- microcosm was characterized by relatively low diversity of Methylococcaceae, with sequences closely related to those of the characterized Methylobacter species being most numerous (Table 3 and Supplemental Table 1). The diversity of Methylophilaceae and Methylocystaceae was also low in this microcosm. Community diversity in the +O₂-NO₃^-microcosm was somewhat higher, and most of the Methylococcaceae sequences were novel sequences that could not be affiliated with any described Methylococcaceae. The Methylophilaceae sequences were also most diverse in this microcosm, including unclassified Methylophilaceae. In the -O₂-NO₃^- microcosm, only two methylotroph sequences were detected, both Methylococcaceae, while both Methylococcaceae and Methylophilaceae were detected in the -O₂+NO₃^- microcosm.

The diversity of methane oxidizing bacteria was further assessed by profiling genes encoding subunits of both particulate (pmo) and soluble (mmo) methane monooxygenase enzymes. Profiling of the pmoB genes (the largest and the less conserved of the pmo genes) revealed significant diversity at the genus level suggesting the presence of at least five identifiable genera within Methylococcaceae and two genera within Methylocystaceae (Table 4), with an exception of the -O₂-NO₃^- microcosm in which no pmoB genes were detected. The Methylocystaceae sequences were only detected in the aerobic microcosms, in agreement with the 16S rRNA gene profiling data. Once again, shifts in relative presence of different phyla were noted suggesting community dynamics in response to variable oxygen tensions and nitrate presence. The lowest phylotype diversity was observed for the -O₂+NO₃^- condition, where only two methanotroph phylotypes were detected, Methylobacter and Methylovulum. In the aerobic microcosms, along with identifiable phylotypes, novel phylotypes of Methylococaceae were present, making up a significant proportion of the population. In the -O₂+NO₃^- microcosm as well as in the unamended sample, pmo sequences were also identified belonging to the NC10 phylum, but in each case these constituted a minor proportion of the population. Very few mmo genes were detected in the metagenomes (Supplemental Table 2), suggesting that most of the active methane oxidizers were devoid of soluble methane monooxygenase.

Genes encoding formaldehyde activating enzymes (fae) were profiled to evaluate the communities possessing formaldehyde oxidation potential. As expected, the diversity of fae genes/proteins extended beyond the bona fide methylotroph species (Table 5) and included Planctomycetes, Archaea, Burkholderiales as well as unclassified bacteria whose ability to oxidize or assimilate C1 compounds is unknown. The highest diversity of Fae was observed in the unamended sample and in the +O₂-NO₃^- microcosm. In the latter, however, the non-methylotroph sequences made up a very minor proportion of total sequences while in the former they made up 30.3% of total sequences. In the enriched microcosms, with the exception of the -O₂+NO₃^- condition, the majority of sequences were classified as belonging to Methylococcaceae and Methylophilaceae. The Methylocystaceae sequences were only prominently present in the +O₂-NO₃^- microcosms. Only in the -O₂+NO₃^- microcosm, the bona fide methylotroph fae genes constituted a relatively minor fraction of the total, with the dominant type being the planctomycete type. While Planctomycetes typically encode Fae and other functions of the tetrahydromethanopterin-linked formaldehyde oxidation pathway, methylotrophy has not been demonstrated in these species, and their genomes lack recognizable genes for methane oxidation or methanol oxidation functions (Chistoserdova, 2011b). The Methylococcaceae were most prominently represented by the Methylobacter and Methylomonas-like sequences, while within Methylophilaceae, the Methylotenera sequences were most prominent. Overall, the community structure predictions as deduced from Fae analysis agreed with other analyses (16S rRNA gene, pmoB and mmoX genes and other methylotrophy genes, not shown).

Methanol is the primary product of methane oxidation, and thus methanol dehydrogenase is an essential enzyme in the methane oxidation pathway (Anthony, 1982). In methanotrophs, this reaction is carried out by an enzyme encoded by mxaFI genes (Chistoserdova & Lidstrom, 2013). mxaFI genes are also essential in methanol oxidation by other groups of methylotrophs (Chistoserdova & Lidstrom, 2013). However, recently methylotrophs capable of methanol oxidation were described not possessing mxaFI genes, and in these, either mdh2 or xoxF genes have been implicated in this function (Chistoserdova & Lidstrom, 2013). Previous metagenomic analysis of Lake Washington populations suggested that the abundant Methylophilaceae phylotypes lacked mxaFI genes but possessed multiple copies of xoxF (Kalyuzhnaya et al., 2008). Genomes of both alpha- and gammaproteobacterial methanotrophs are also known to encode xoxF genes (Chistoserdova, 2011b). The metagenomes generated in this study were specifically analyzed for the presence of mxaF and xoxF genes that classed in the families of Methylococcaceae, Methylocystaceae and Methylophilaceae, using high stringency BLAST searches (Supplemental Table 3). In each category, both types of genes were detected. As previously observed (Kalyuzhnaya et al., 2008), the Methylophilaceae populations appeared to be dominated by the types lacking mxaF. The Methyococcaceae xoxF genes also appeared to be more abundant compared to the mxaF genes with the exception of the +O₂-NO₃^- microcosm. This may suggest either that some of the Methylococcaceae lack mxaF genes or that some of them contain multiple copies of xoxF genes.

The diversity of genes for denitrification and relative abundance of genes belonging to key methylotroph groups (the Methylococcaceae, the Methylophilaceae and the Methylocystaceae) was evaluated by taxonomic profiling of the genes identified by the IMG/M platform as nitrate reductase, nitrite reductase, nitric oxide reductase or nitrous oxide reductase genes. For comparison, diversity and relative abundance of methylotroph nitrogenase genes were evaluated, revealing a significant phylogenetic complexity of the communities with a potential for denitrification (Table 6 and Supplemental Tables 4–8). We found that in the aerobic microcosms, the nitrogenase, nitrate reductase and nitrite reductase genes taxonomically ascribed to Methylococcaceae and primarily to the genus Methylobacter, were most abundant (Fig. 4), suggesting that Methylobacter was the major species capable of both nitrogen fixation and denitrification in these conditions, as well as in the native sediment. Methylophilaceae-affiliated nitrate- and nitrite reductases represented another dominant class in the aerobic microcosms (Fig. 4). The most abundant nitric oxide reductase type in these microcosms was classified as Methylophilaceae and most prominently Methylotenera, suggesting that most of Methylotenera species encode this step of the pathway while some of the Methylococcaceae represented in the metagenomes lack the respective gene. The phylotypes classified as Methylobacter were notably under-represented in this category (Supplemental Table 6). No nitrous-oxide reductase genes were detected affiliated with either Methylococcaceae or Methylophilaceae suggesting that the denitrification pathway may be incomplete in these species. However, a number of sequences were classed with Methylocystaceae in the +O₂-NO₃^- microcosm (Fig. 4). Analysis of the nitrogenase genes demonstrated that Methylococcaceae and most prominently Methylobacter species were the major type possessing the potential for nitrogen fixation (up to 80% of total sequences), with a prominent presence of Methylocystaceae in the +O₂-NO₃^- microcosm (Fig. 4, Supplemental Table 8). The remaining sequences of the denitrification and nitrogen fixation genes were distributed evenly among a variety of phyla, and no other dominant groups or groups specifically responding to nitrate were detected (Supplemental Tables 4–8).