General description of sequenced metagenomes and bioinformatic strategy
A metagenomic approach was used to characterize different anaerobic communities, and to assess their ability for methane production. Among the native communities, two originated from full-scale biogas reactors. Samples were collected from an agricultural biogas fermenter (ABF) and an agricultural biogas hydrolyzer (ABH) and were treated as a reference source of an efficient methanogenic community residing in a two-stage biogas plant. Aside from the extensively studied biogas reactor communities, we examined also more natural communities from environments where methane emission is detected. These other communities came from cattle manure (CM), cattle slurry (CS) and sewage sludge from a wastewater treatment plant (WTP), as biogas plants are often supplemented with these waste materials. The last two environmental samples originated from a lowland bog (LB) and an ancient gold mine (GM), as these habitats are a reservoir of diverse uncultivated microorganisms potentially beneficial to the methanogenesis process [21,22,23]. For comparative analyses of microbial community structural and functional change following a cultivation process, two laboratory consortia (ABF_TS and WTP_TS) were also implemented in the study. They were selected in batch reactors, and then stabilized in two-stage (TS) bioreactors.
Altogether, nine samples were analyzed in this study. Basic physico-chemical parameters were determined for all native and laboratory communities (see Additional file 1: Table S1) and metagenomic DNA was isolated and deep sequenced on the Illumina HiSeq 1500 platform. We obtained approximately 10–13 gigabases (79 million sequence counts on average) for each sample. The statistics of sequence counts and annotation analysis are presented in Additional file 1: Table S2. The reads obtained were analyzed using two different approaches implemented in the most commonly used pipelines (Fig. 2).
General analyses of community taxonomic and functional profiles were carried out through identification of similar sequences in reference databases with the use of the MG-RAST pipeline [24]. To obtain a more detailed insight into the methanogenesis process and the microbes involved, we used MetAnnotate [25], in which a HMM-based search is implemented and taxonomy is determined by the best hit or phylogenetic tree placement of obtained hits. In this study, we compared results from both approaches.
Taxonomic profiling
Taxonomic diversity was determined from metagenomic sequences annotated against the RefSeq database using the MG-RAST pipeline [24] (see Additional file 1: Table S3 and Tax MG-RAST, Additional file 2). At the domain level, all sequenced microbial communities were dominated by Bacteria (87.4–98.9%), followed by Archaea (0.4–11.5%) and Eukaryota (0.5–1.3%). Community structure analysis at the genus level showed 970 to 1525 different taxonomic groups in the samples. The most diverse and even was community originated from LB sample (Shannon–Wiener index—5.795, Pielou index—0.811) whereas the least versatile and uniform was ABH community (Shannon–Wiener index—4.055, Pielou index—0.590) (Additional file 1: Table S4). However in all samples there were only two to eight genera with the relative abundance above 2% (Fig. 3a). Microbes of Bacteroides and Clostridium genera were among the most common and the most abundant (Fig. 3b).
Among the investigated metagenomes, samples coming from agricultural biogas plant (ABF, ABH) were considered model consortia that harbor microorganisms essential for all steps of high-performance biogas production. The most specific feature of the agricultural biogas plant fermenter sample (ABF) was the high abundance of methanogenic archaeons, i.e. Methanoculleus (4.8%) and Methanosarcina (2.3%), whereas in other native metagenomes their abundance was significantly lower, close to 0.7% (ABH and LB samples) or even not exceeding 0.1% (CM, CS, GM and WTP samples (see Additional file 1: Table S3). Furthermore, in the agricultural biogas plant samples, we detected 2.6% (ABF) and 3.5% (ABH) of metagenomic sequences classified as Candidatus Cloacamonas, while in other native metagenomes sequence counts of this bacterium accounted for a maximum of 0.2%. Moreover, the ABH sample was highly enriched in polysaccharide-degrading Prevotella (28.5%), Lactobacillus (4.1%) and numerous genera belonging to Veillonellaceae family (2.2% cumulatively).
Interestingly, in the cattle-associated samples we observed that Pseudomonas (7.5%; 6.7%; CM and CS, respectively), Acinetobacter (1.5%; 4.2%) and Bacillus (2.4%; 2.2%) were more abundant than in bioreactor samples. In addition, CM was richer in members of Acholeplasma (2.4%). In the case of the community originated from a wastewater treatment plant (WTP), we observed a high number of sequences assigned to Acinetobacter (9.6%), Arcobacter (6.1%), Aeromonas (4.6%) and Acidovorax (3.3%). The last two environmental samples analyzed, originating from a gold mine (GM) and a lowland bog (LB), were the most diverse and had few common genera with each other and with the remaining metagenomes (Fig. 3; Additional file 1: Fig. S1). The most unique feature of GM sample was a very high content of sequences assigned to methanotrophic Methylobacter (11.3%), Methylococcus (2.5%) and Methylotenera (2.3%). As a comparison, in other metagenomes these bacteria accounted for a maximum of 0.5%. In the case of LB community, the highest abundance was observed for Geobacter (3.3%) and Candidatus Solibacter (2.4%) (Fig. 3b). Importantly, LB community was enriched in numerous, low abundance environmental Archaea, including both methanogenic Archaea (3.9%) and other Archaea (4.7%) without a confirmed methane production ability (see Tax MG-RAST, Additional file 2).
Both laboratory consortia, originated from an agricultural biogas fermenter (ABF_TS) and sewage sludge from a wastewater treatment plant (WTP_TS), were dominated by Bacteria (93.0% in ABF_TS; 94.6% in WTP_TS), followed by Archaea (6.1% in ABF_TS, 4.9% in WTP_TS). Taxonomic annotations of the ABF_TS sample revealed that despite laboratory cultivation, the structure of microorganisms has been largely preserved (Fig. 3b; Additional file 1: Fig. S1). Nevertheless, we observed at least 1% change of the relative abundance of the following dominant genera: Bacteroides (from 10.9 to 14.1%), Clostridium (from 6.0 to 8.5%) and Parabacteroides (from 2.2 to 3.3%), Pseudomonas (from 0.3 to 1.3%), Prevotella (from 4.4 to 2.7%) and Methanoculleus (from 4.8 to 0.5%) (Fig. 3b; Additional file 1: Table S3). In the case of the consortium selected from sewage sludge of a wastewater treatment plant (WTP_TS), a substantial change of the community structure was observed compared to the native WTP sample (Fig. 3b; Additional file 1: Fig. S1). In general, abundance of Proteobacteria representatives diminished from 56.7 to 11.7%, while abundance of Bacteroidetes and Firmicutes increased from 26.8 and 10.0% in inocula to 34.6 and 33.8% in the laboratory consortium, respectively. Moreover, Archaea abundance increased from 0.5 to 4.9% of the total microbial structure. It seems that the most important change of laboratory consortium in comparison to its inoculum was the increase of fermentative Bacteroides (from 11.5 to 19.5%), Clostridium (from 2.1 to 10.7%), Parabacteroides (from 2.0 to 3.3%), Eubacterium (from 0.8 to 2.1%), Candidatus Cloacamonas (from 0.1 to 3.1%), Ruminococcus (from 0.4 to 1.6%), Syntrophomonas (from 0.1 to 1.6%) and methanogenic Methanosarcina (from 0.1% to 3.1%) (see Additional file 1: Table S3). Simultaneously abundance of other genera such as Acinetobacter, Arcobacter, Aeromonas, Acidovorax, Streptococcus, Flavobacterium, Thauera and Tolumonas decreased from 9.6% in native community (WTP) to max. level of 0.6% in laboratory community (WTP_TS) (Fig. 3b; Additional file 1: Table S3).
Functional profiling
To explore the metabolic potential of the studied communities, we performed a detailed analysis of metagenomic sequences annotated against SEED subsystems within the MG-RAST pipeline. We detected on average 5590 functional categories with at least 0.001% annotated reads (see Additional file 1: Table S2 and Fun MG-RAST, Additional file 2). Analysis of commonly used diversity and evenness indices showed that all metagenomes were quite similar in terms of diversity. The most functionally diverse and even was WTP community (Shannon–Wiener index—7.560, Pielou index 0.834) while the least one was ABF metagenome (Shannon–Wiener index—7.196, Pielou index 0.816) (see Additional file 1: Table S5). Furthermore, similarly to RefSeq Bray–Curtis distances calculation, some samples located near each other e.g. ABF and ABF_TS, while GM and LB samples were the most different from majority of the analyzed metagenomes (see Additional file 1: Fig. S2).
We then focused on functional analysis relevant to methanogenesis pathways. Cumulative relative abundances of genes involved in various methanogenesis pathways were in the range of 0.31–1.15% (Fig. 4), with a caveat that genes from acetoclastic and methylotrophic pathways can contribute to processes other than methanogenesis.
Considering the percentage of methanogenesis-related annotations, ABF community followed by LB, ABF_TS and WTP_TS and ABH consortia had the highest potential for methane production via different pathways. Communities from an agricultural biogas plant and laboratory reactors (ABF, ABH, ABF_TS, WTP_TS) had similar sequence profiles of the selected genes. However, the functional profile of lowland bog (LB) metagenome differed significantly (Fig. 4). In the LB sample, we observed an overrepresentation of acs, cdh and hdr genes compared to other methanogenesis genes, including the key mcr genes. An even greater overrepresentation of genes involved in acetate utilization was observed for cattle-associated (CM, CS) and wastewater treatment plant (WTP) samples. However, in these samples most of the selected sequences were classified to ack and pta genes but few to cdh genes that directly mediate the transfer of carbon to methane pathway (Fig. 4). In the case of the native consortium from a gold mine (GM), we observed a relatively high number of sequences of fmd, ftr, mch, hmd, ack, pta, acs genes involved in the utilization of carbon dioxide and acetate, and very few sequences of mcr genes encoding the key methyl-CoM reductase enzyme.
The comparison of laboratory consortia and their native counterparts indicated that during laboratory cultivation the cumulative abundance of methanogenesis genes increased for WTP_TS and decreased for ABF_TS. In the case of WTP_TS, the increase concerned all genes except ack, pta and acs, while for ABF_TS the reduction concerned all studied genes.
Methanogenesis genes-specific phylogenetic characterization
A sequence-based analysis of metagenomic sequences with the use of the MG-RAST pipeline gives an opportunity to explore metabolic potential of complex communities. Nevertheless, metagenomic data enables also a simultaneous identification of function and microorganisms responsible for specific processes. In the following part of our study, we present detailed results of phylogenetic placement of methanogenesis-related sequences. For this type of an analysis, the HMM search and the taxonomic classification approach implemented in MetAnnotate [25] were used. We focused on genes encoding methyl CoM reductase (mcr) responsible for the final release of methane, as well as on sequences of genes encoding enzymes which are considered crucial in the utilization of various substrates, such as: CO2—formylmethanofuran dehydrogenase (fmd); CH3COOH—CO dehydrogenase/acetyl-CoA synthase (cdh); CH3OH—methanol methyltransferase (mta); (CH3)xNH—methylamine methyltransferase (mtm, mtb, mtt). A summary of key microorganisms involved in specific steps of the methanogenesis superpathway for each analyzed sample is presented in Additional file 1: Fig. S3A–I.
The Metannotate-based analysis of native communities from an agricultural biogas plant (ABF and ABH) revealed that they were highly enriched in methyl-CoM reductase (mcr) genes. A high abundance of mcr genes was also observed in laboratory consortia (ABF_TS and WTP_TS) compared to the other environmental communities, such as CM, CS, GM, LB, WTP (Fig. 5a, bar graph panel). This is consistent with the data obtained using MG-RAST (Fig. 4). Metannotate-based analysis of mcr genes revealed substantial differences between analyzed samples in their phylogenetic placement. As was showed on Fig. 5a, most of the mcr reads of ABF sample were assigned to hydrogenotrophic Methanoculleus (59%), while in ABH they were dominated by Methanobacterium (39%) and Methanoculleus (30%). The largest number of mcr sequences in both of the laboratory communities (ABF_TS and WTP_TS) mapped to Methanosarcina (17%, 47%), Methanoculleus (29%, 16%) and Methanocorpusculum (22%, 16%), respectively. In the remaining metagenomes, most of the mcr sequences mapped also to hydrogenotrophic methanogens. In the cattle-related samples, Methanobrevibacter (23% in CM, 38% in CS) and Methanocorpusculum (31% in CM, 13% in CS) dominated, whereas Methanoculleus (20%) and Methanoregula (22%) dominated in GM and LB samples, respectively (Fig. 5a). In contrast, in WTP consortium, we observed a high number of mcr sequences mapped to acetoclastic Methanosaeta (17%). We also detected a substantial number of hits to representatives of the seventh order of methanogens, such as Candidatus Methanomethylophilus (24% in CM, 22% in CS) and Methanomassiliicoccus (12% in LB, 12% in WTP) (Fig. 5a).
Apart from the identification of microorganisms responsible for the last step of methanogenesis, we aimed for identifying those involved in carbon incorporation from different substrates, such as CO2, CH3COOH, CH3OH and (CH3)xNH. An analysis of formylmethanofuran dehydrogenase subunits ACE (fmdACE), indicators of hydrogenotrophic pathway, showed that among the native communities most hits were detected for ABF consortium, followed by GM, LB, ABH, ABF_TS, WTP_TS, CM, CS and WTP communities (Fig. 5b, bar graph panel). There was an up to 19-fold difference in the fmdACE abundance between the ABF and WTP samples. In the ABF sample, hydrogenotrophic Methanoculleus (66%) was the organism predominantly assigned to fmd genes. In the WTP sample, a majority of fmd sequences mapped to Thauera (21%) and to numerous taxa of Archaea and Bacteria, which did not exceed 6% (Fig. 5b). Following the laboratory cultivation of the aforementioned communities, Methanosarcina predominated in ABF_TS (28%) and WTP_TS (54%) fmd sequences. In comparison, for LB metagenome many of the fmd sequences mapped to various low-abundant Archaea and Bacteria with the highest counts for Desulfobacter (7%) and Methylobacterium (7%). Interestingly, in the case of GM metagenome, most fmd sequences were classified to methanotrophic bacteria, such as Methylobacter, Methyloglobulus, Methylovulum (12% each), but few were classified to methanogenic Archaea (Fig. 5b). In contrast, agricultural biogas hydrolyzer (ABH), cattle manure (CM) and cattle slurry (CS) were dominated mainly by hydrogenotrophic methanogens, such as Methanoculleus (35%) and Methanobacterium (32%) in ABH sample, Methanocorpusculum (23%) in CM sample and Methanobrevibacter (22%) in CS sample. Furthermore, both CM and CS had many hits (~ 10%) for Methanospirillum, Methanoregula, Anaerosalibacter and Blautia (Fig. 5b).
In order to identify microorganisms involved in the acetoclastic pathway, we analyzed phylogenetic assignments of the D subunit of CO dehydrogenase/acetyl-CoA synthase (cdhD), as it is directly involved in the transmission of a methyl group from acetate during acetoclastic methanogenesis [4, 5]. The abundance of cdhD was highest for ABH, followed by LB, ABF, ABF_TS, WTP_TS, CM, CS samples, while in GM and WTP metagenomes it was scarce (Fig. 5c, bar graph panel). Identification of phylogenetic matches of cdhD indicated that in most of the metagenomes acetate utilization was mediated by two to four dominant genera with percentage above 10%. In the case of laboratory communities, cdhD sequences were mostly assigned to one methanogen, namely Methanosarcina, which accounted for 65% (ABF_TS) and 94% (WTP_TS). In the case of the other two samples, approximately half of the annotations were classified to one methanogen, namely Methanobacterium (53%) in ABH and Methanoregula (49%) in LB (Fig. 5c). For other metagenomes, cdhD sequences were assigned to few dominant genera. These were Methanosarcina (32%), Methanoculleus (24%), Methanosaeta (17%), and Methanobacterium (12%) in the ABF metagenome; Methanoregula (24%) and Desulfocapsa (15%) in GM; Methanosaeta (31%) and Methanoregula (23%) in WTP and Methanosarcina (24%, 32%), Blautia (25%, 15%), Treponema (13%, 14%) in CM and CS samples, respectively (Fig. 5c). Furthermore, in the CM, CS, GM, LB, WTP samples, cdhD was often encoded by various Archaea and Bacteria.
Utilization of methanol or methylamines is the third commonly recognized methanogenic pathway, which contains genes of methanol and mono-, di- and trimethylamine methyltransferases (mta, mtm, mtb, mtt, respectively). Comparison of the available domain profile sequences of mtaB, mtmB, mtbB and mttB showed that their abundances were almost constant between metagenomes, with the highest count for laboratory WTP_TS community and the lowest count for the native WTP consortium. The GM and LB metagenomes were an exception, as for them we detected enrichment of mttB compared to the other known methylotrophic domain profiles. Phylogenetic matches of functional annotations for genes of methylotrophic pathway showed that the majority of agricultural biogas plant, cattle-associated and laboratory sample sequences were assigned to Archaea, such as Methanosarcina, Candidatus Methanomethylophilus and Candidatus Methanoplasma (Fig. 6). The predominance of Methanosarcina (up to 99%) was particularly apparent for agricultural biogas plant and laboratory samples (ABF, ABH, ABF_TS, WTP_TS), while representatives of Methanomassiliicoccales were more abundant in CM and CS, for example, up to 60% of sequences of Candidatus Methanomethylophilus. In the case of GM, LB and WTP communities, we observed a higher contribution of bacterial genera, such as Diplosphaera, Desulfococcus, Levilinea and Eubacterium (even up to 66% of all annotations) (Fig. 6).
Amplicon-based analysis of methanogenesis-related genes
In order to test if methanogenesis potential of deeply sequenced communities could be grasped by widely used and simple amplicon-based approach, we performed a semi-quantitative amplification of marker genes designed to monitor the methanogenesis potential of environmental samples [18]. PCR-amplified genes included alpha, beta and gamma subunits of methyl-CoM reductase (mcrABG), methanol-specific methyltransferase (mtaB) and methylamine-specific methyltransferase (mtbA).
Results showed that all five analyzed genes were successfully amplified from total DNA originating from ABF, ABH, CS, CM, WTP and laboratory communities ABF_TS and WTP_TS, which suggested their potential to carry out methane fermentation. However, lowland bog (LB) and gold mine (GM) communities seemed to have weak methanogenesis potential, as for the LB sample only mcrA product and a low amount (a faint band) of mtaB product were obtained and for the GM sample none of the PCR products were obtained (Fig. 7). In contrast, the control amplification of bacterial and archaeal 16S rDNA fragments was positive for all samples (Additional file 1: Fig. S4).