Metabolite measurements

Concentrations of acetate, propionate, and butyrate in the liquid sample of the soil slurries were measured by HPLC equipped with an ion-exclusion column (Aminex HPX-87-H, BioRad, München, Germany) and coupled to a UV–Vis detector (Sykam, Fürstenfeldbruck, Germany) [38]. In addition, gas samples were taken from the same set of slurries for process measurements. A GC-8A gas chromatograph (Shimadzu, Duisburg, Germany) containing a Haysep Q column was used to measure CH₄ and CO₂. Data were analyzed with PeakSimple software (SRI Instruments, Bad Honnef, Germany) and calculated by linear regression [12, 36, 39].

DNA and RNA extraction

Total DNA was extracted from soil slurries using the DNeasy® PowerLyzer® PowerSoil kit (QIAGEN, Hilden, Germany) according to the manufacturer’s instructions. Agarose gel (1%) electrophoresis and fluorometry were used to check for the integrity and quantity of each DNA extract. Fluorometric measurements were done on a Qubit 2.0 Fluorometer using the Qubit dsDNA BR Assay Kit (Thermo Fisher Scientific, MA, USA).

The RNeasy® PowerSoil Total RNA Kit (QIAGEN, Hilden, Germany) was used for the extraction of total RNA. The extraction procedure followed the manufacturer’s instructions. The RNA extracts were treated with DNase I (Ambion, Austin, USA) and purified using the RNA Clean and Concentrator kit (Zymo Research, CA, USA). The integrity of purified RNA was assessed using the Experion^TM RNA HighSens Analysis Kit (Bio-Rad, CA, USA), while the yield was determined with the Qubit^TM RNA HS Assay Kit (Thermo Fisher Scientific, MA, USA).

Quantitative PCR (qPCR) and reverse transcription qPCR (RT-qRCR)

Primer sets and temperature profiles for quantitation of bacterial 16S rRNA and mcrA genes (qPCR) and their transcripts (RT-qPCR) are shown in Table S1. The mcrA gene is a standard biomarker to detect and quantify methanogens in environmental samples [36, 39, 40]. In RT-qPCR, randomly reverse-transcribed RNA was generated using the GoScript Reverse Transcription System (Promega, Mannheim, Germany) according to the manusfacturer’s instructions. The standard curve for quantifying bacterial 16S rRNA genes and transcripts was constructed using the genomic DNA of Escherichia coli (calibration range from 10 to 10⁹ copies). The standard curve for quantifying mcrA genes and transcripts was constructed using a mcrA fragment cloned into the pGEM®-T Easy plasmid (Promega, WI, USA) (calibration range from 10 to 10⁸ copies) [39]. The mcrA fragment used for cloning was obtained from the genomic DNA of Methanosarcina barkeri [41]. The qPCR reactions were carried out on an iCycler Real-Time PCR Detection System (CFX Connect^TM, Bio-Rad). The PCR efficiency was at least 85% (R² > 0.98). The presence of unspecific products was checked by melt curve analysis.

Illumina library preparation and sequencing

Libraries for sequencing were prepared for both metatranscriptomic (cDNA) and metagenomic (total DNA) analysis.

Computational analysis of metatranscriptomic datasets

Computational analysis of metagenomic datasets

qPCR and RT-qPCR

Genes and transcripts of bacterial 16S rRNA and methanogenic mcrA were quantified over the 120-day incubation period (Fig. ​(Fig.2,2, Tables S6, S7, S8, S9). Changes in the copy numbers of bacterial 16S rRNA genes and transcripts (Fig. ​(Fig.2a,2a, c) and mcrA genes and transcripts (Fig. ​(Fig.2b,2b, d) showed an M-like up and down over incubation time. The copy numbers of bacterial 16S rRNA genes and transcripts peaked first on day 21 (1.2 × 10¹⁰ genes vs. 4.8 × 10¹¹ transcripts g^-1 dry soil) and again on day 35 (1.2 × 10¹⁰ genes vs. 4.2 × 10¹¹ transcripts g^-1 dry soil), with an intermediate decrease on day 28 (5.7 × 10⁹ 16S rRNA genes vs. 2.7 × 10¹¹ transcripts g^-1 dry soil) (Fig. ​(Fig.22a,c).

Open in a separate window

Fig. 2

Copy number quantification of bacterial 16S rRNA and methanogenic mcrA genes (qPCR) and their transcripts (RT-qPCR) in Philippine paddy soil slurries over the 120-day incubation period. Copy numbers of bacterial 16S rRNA genes (a) and transcripts (c) per gram of dry paddy soil. Copy numbers of mcrA genes (b) and transcripts (d) per gram of dry paddy soil. Data are means ± SE (n = 3). Differences between two incubation time points were determined using one-way ANOVA (P_FDR<0.05). The exact copy numbers per gram of dry paddy soil are shown in Tables S6, S7 (genes) and S8, S9 (transcripts)

Likewise, the copy numbers of mcrA genes and transcripts peaked first on day 21 (2.9 × 10⁹ mcrA genes vs. 9.2 × 10¹⁰ transcripts g^-1 dry soil) and again on day 35 (3.4 × 10⁹ mcrA genes vs. 6.5 × 10¹⁰ transcripts g^-1 dry soil), with an intermediate decrease on day 28 (1.3 × 10⁹ mcrA genes vs. 3.8 × 10¹⁰ transcripts g^-1 dry soil) (Fig. ​(Fig.2b,2b, d).

After day 35, the gene and transcript copy numbers of both bacterial 16S rRNA and mcrA strongly decreased towards day 120 (Fig. ​(Fig.22).

The bacterial metatranscriptome (16S rRNA, mRNA)

Taxon-specific dynamics

The family-level population dynamics varied over incubation time (Figs. ​(Figs.33 and S3, Tables S10 and S11). At the rRNA level, the relative abundance of Geobacteraceae peaked between days 7 and 11 (Fig. ​(Fig.3a).3a). Thereafter, the abundance of Clostridiaceae and Peptococcaceae peaked at day 21, while Lachnospiraceae showed the maximum abundance on day 28 (Fig. ​(Fig.3a,c).3a,c). In the late phase (days 35 to 120), Heliobacteriaceae and Anaerolineaceae showed increased rRNA abundances and were, together with Geobacteraceae, the predominant bacterial populations (Fig. ​(Fig.33a,c).

Open in a separate window

Fig. 3

Relative abundance changes of dominant bacterial families (> 2%) on rRNA (a, c) and mRNA (b, d) levels, respectively. The percentage abundances are given in relation to total bacterial 16S rRNA and mRNA, respectively. Data are means ± SE (n = 3). Asterisks * (P_FDR≤0.05), ** (P_FDR≤0.01), and *** (P_FDR≤0.001) indicate significant differences. Asterisks shown aside from the taxonomic names indicate that in STAMP, the relative rRNA or mRNA abundances significantly changed through the whole incubation period. Asterisks directly shown in the plots indicate that in DESeq2, the rRNA or mRNA abundance in that particular sampling time point significantly differed from the original abundance at the first sampling (day 3) or vice versa that the relative rRNA or mRNA abundance at day 3 (Symbiobacteriaceae) significantly differed from all other sampling time points. Significance values (P_FDR) are only indicated on rRNA and mRNA levels for the taxon-specific abundance peak. The complete set of P_FDR values obtained by DESeq2 analysis for all sampling time points is shown in Tables S10 (rRNA) and S11 (mRNA)

Like on the 16S rRNA level, the mRNA abundance of the Geobacteraceae peaked on day 7, covering 22% of total community-wide mRNA and 36% of total bacterial mRNA. This was followed by a rapid abundance decline towards day 21 and another increase thereafter (Figs. ​(Figs.3b and3b and S4). The Peptococcaceae reached their peak mRNA abundance on day 21 (Fig. 3d). Likewise, Clostridiaceae (day 21) and Lachnospiraceae (day 28) showed peak mRNA abundances consistent with those observed on the rRNA level. The Acidobacteriaceae and Anaerolineaceae displayed higher mRNA abundances in the late phase (days 35 to 120) than in the early phase (days 3 to 21) (Fig. ​(Fig.33d). In addition, putative syntrophic populations, such as Peptococcaceae and Heliobacteriaceae, reached a second peak abundance during the late phase (Fig. ​(Fig.3b,3b, d).

The methanogen metatranscriptome (16S rRNA, mRNA)

Methanosarcinaceae, Methanocellaceae, and Methanotrichaceae were the prevailing methanogenic families throughout the complete incubation period (Fig. ​(Fig.4a,4a, b; Tables S10 and S11). The Methanosarcinaceae was the dominant methanogen group, showing two-peak abundance dynamics. The family-level rRNA and mRNA abundances peaked first around days 11 and 14, then decreased until day 28, but increased thereafter again and reached a second activity peak between days 35 and 60. The Methanocellaceae were first detectable on day 11. Following the peak abundance on day 21, the relative rRNA and mRNA abundances steadily decreased towards day 120. Significant transcript levels of the Methanotrichaceae were only detectable on mRNA level during the late phase, with the peak abundance at day 120.

Open in a separate window

Fig. 4

Relative abundance changes of the dominant methanogenic families (> 2%) on rRNA (a) and mRNA (b) levels, respectively. Their percentage abundances are given in relation to total archaeal 16S rRNA and mRNA, respectively. Data are means ± SE (n = 3). Asterisks * (P_FDR≤0.05), ** (P_FDR≤0.01), and *** (P_FDR≤0.001) indicate significant differences. Asterisks shown aside from the taxonomic names indicate that in STAMP, the relative rRNA or mRNA abundances significantly changed through the whole incubation period. Asterisks directly shown in the plots indicate that in DESeq2, the rRNA or mRNA abundance in that particular sampling time point significantly differed from the original abundance at the first sampling (day 3). Significance values (P_FDR) are only indicated on rRNA and mRNA levels for the taxon-specific abundance peak. The complete set of P_FDR values obtained by DESeq2 analysis for all sampling time points is shown in Tables S10 (rRNA) and S11 (mRNA)

Defining the dominant Methanosarcina populations

The assembled near-full-length 16S rRNA sequences grouped into four distinct Methanosarcina populations (Groups I to IV) (Fig. ​(Fig.5a).5a). These differed in their abundance dynamics (Fig. ​(Fig.5b).5b). The Group II population was most closely related to Methanosarcina sp. MSH10X1 and predominant throughout slurry incubation. Methanosarcina Group IV was the second most abundant population and closely affiliated with M. horanabensis. The Group II and Group IV populations showed opposite 16S rRNA abundance dynamics, with both Group II minimum abundance and Group IV maximum abundance at day 11 (51% [Group II] versus 34% [Group IV] of total Methanosarcina 16S rRNA) (Fig. ​(Fig.5b).5b). The Methanosarcina Group I and Group III populations displayed low but relatively stable transcript abundances (collectively < 20% of total Methanosarcina 16S rRNA) throughout slurry incubation (Fig. ​(Fig.55b).

Open in a separate window

Fig. 5

Phylogenetic tree of the four Methanosarcina populations detected in the paddy soil slurries (a) and their relative abundance dynamics over incubation time (b). a The neighbor-joining tree was constructed based on near full-length 16S rRNA sequences (> 1200 bp) assembled by EMIRGE from the metatranscriptomic datasets (415 sequences) and reference sequences extracted from Methanosarcina genomes in GTDB. b The relative abundance dynamics of the four distinct Methanosarcina populations (Groups I to IV) were inferred from the assembled 16S rRNA sequences. Data are means ± SE (n = 3)

Mapping-independent expression analysis of methanogenic pathways

The collective mRNA abundance of major methanogenesis pathways (Fig. ​(Fig.6a6a and Table S13) and the transcript abundance of particular pathway marker genes (Fig. ​(Fig.6b6b and Tables S13, S14) varied with incubation time. Acetoclastic methanogenesis was the dominant methane production pathway with the greatest transcript abundance in the early phase (36.4% of total mRNA assigned to KEGG level 3 category “methane metabolism”). This involved the peak transcript abundance of the following pathway marker genes: cdhAB, cdhCDE, ack, and pta (Fig. ​(Fig.6b).6b). The changes in their relative expression level agreed well with the overall mRNA abundance dynamics of the Methanosarcinaceae during the early phase (Figs. ​(Figs.4b,4b, ​b,6 a-c).6 a-c). Transcripts of genes (acs) indicative of Methanotrichaceae-driven acetoclastic methanogenesis steadily increased in relative abundance during the late phase. Their peak abundance occurred on day 120, at which the relative expression level of genes indicative of Methanosarcinaceae-driven acetoclastic methanogenesis had significantly declined (Figs. ​(Figs.4b4b and 6b, c).

Open in a separate window

Fig. 6

Relative mRNA abundance dynamics of individual methanogenic pathways over incubation time a, transcript dynamics of related key pathway genes b, and taxonomic assignment of the transcripts c. A list of the full gene names is shown in Table S14. The relative abundance values are given in relation to total mRNA affiliated to the KEGG level 3 category “methane metabolism”. The relative expression levels were calculated based on TPM values. Data are means ± SE (n = 3). Asterisks * (P_FDR≤0.05), ** (P_FDR≤0.01), and *** (P_FDR≤0.001) indicate significant differences. Asterisks shown aside from the gene names indicate that in STAMP, the relative mRNA abundance significantly changed through the whole incubation period. Asterisks directly shown in the heatmap (b) indicate that in DESeq2, the mRNA abundance in that particular sampling time point significantly differed from the original abundance at the first sampling (day 3) or vice versa that the relative mRNA abundance (ack) at day 3 significantly differed from all other sampling time points. Significance values (P_FDR) are only indicated for the peak transcript abundance(s) of each key pathway gene. c The percentage of acs transcripts affiliated with Methanocellaceae (acetate assimilation) was 68.1% on day 120, while those affiliated with Methanotrichaceae (acetoclastic methanogenesis) contributed 31.9%. The complete set of P_FDR values obtained in DESeq2 analysis for all sampling time points is shown in Table S13

Members of the Methanocellaceae were the prevailing H₂/CO₂-utilizing methanogens, with genes encoding hydrogenotrophic methanogenesis (fwd, ftr, mch, mtd, mer) being most expressed by this family-level group (Fig. ​(Fig.6b,6b, c). Their gene expression dynamics agreed well with the overall mRNA abundance dynamics of the Methanocellaceae (Figs. ​(Figs.4b and4b and 6a, b).

In addition, transcripts of genes (mtaA, mtaBC, mtbA, mtbBC, mtmBC, mttBC) indicative of methylotrophic methanogenesis were detected, with mtaBC exhibiting the greatest transcript abundance throughout slurry incubation. The relative abundance of transcripts involved in methylotrophic methanogenesis increased from day 3 onwards and reached the first peak abundance (5.8% of total mRNA assigned to KEGG level 3 category “methane metabolism”) on day 14, but then decreased towards day 21 (2.1%). Thereafter, the relative transcript level of methylotrophic methanogenesis increased again and reached a second peak (3.4% of total mRNA assigned to KEGG level 3 category “methane metabolism”) around day 60 (Fig. ​(Fig.6a, b).6a, b). This two-peak (14 and 60 days) abundance pattern of mRNA encoding methylotrophic methanogenesis agreed well with the overall abundance dynamics of Methanosarcina-affiliated rRNA (Fig. ​(Fig.4a)4a) and mRNA (Fig. ​(Fig.4b),4b), in both early and late phases. Indeed, the taxonomic assignment of the transcripts involved in methylotrophic methanogenesis showed that they were entirely expressed by the Methanosarcinaceae (Fig. ​(Fig.66c).

Gene expression of methanogenic pathways by Methanosarcina Group II

Mapping-dependent gene expression analysis

KEGG-annotated transcripts were mapped onto the three MAGs (Figs. ​(Figs.7a7a and S11). On KEGG level 2, transcripts involved in energy metabolism showed the greatest mapping efficiency (54% to 86% of total mapped mRNA) throughout slurry incubation, while transcripts encoding the translational apparatus were mapped with high frequency (30–32% of total mapped mRNA), particularly during the first week of slurry incubation (Fig. S11a and Table S17). On KEGG level 3, transcripts involved in methane metabolism were dominantly mapped (51 to 80% of total mapped mRNA) onto the MAGs throughout slurry incubation (Fig. S11b and Table S18). On KEGG level 4, transcripts involved in acetoclastic methanogenesis and methylotrophic methanogenesis were mapped with high frequency (Fig. ​(Fig.7a7a and Table S19). Overall, the transcript mapping frequency of the methylotrophic pathway (5.1–7.2% of total mapped mRNA) was lower than that of the acetoclastic pathway (7.5–12.4% of total mapped mRNA). Among total methanogen mRNA, mtaABC transcripts showed the greatest mapping frequency (Fig. ​(Fig.77a).

Open in a separate window

Fig. 7

a Relative abundance of pathway-specific mRNA and transcripts of key pathway genes in relation to total mRNA that could be mapped onto the Methanosarcina MAGs. A list of the full gene names is shown in Table S14. Mapping onto the MAGs was conducted for each mRNA dataset retrieved in triplicate from the same incubation time point (days 21, 28, and 35). The small dots indicate the mapping efficiency of each mRNA replicate dataset, while the large dot represents the mapping efficiency averaged across all three mRNA replicate datasets of a given sampling time point. The complete set of P_FDR values obtained in DESeq2 analysis for all sampling time points is shown in Table S19. b Relative transcript dynamics of acetoclastic and methylotrophic methanogen pathways inferred from the competitive transcript mapping across the mix of Methanosarcina Group I to IV reference genomes. The relative mapping efficiencies were calculated based on TPM values. Data are means ± SE (n = 3). The mapping statistics are shown in Table S20

Linking metabolite turnover with methanogen dynamics

Typical intermediates (i.e., acetate, propionate, butyrate) of the methanogenic food web transiently accumulated and were then rapidly consumed. Acetate is the most abundant intermediate during the anaerobic organic matter breakdown [19, 70]. Generally, it is the direct substrate for acetoclastic methanogens, including Methanosarcinaceae and Methanotrichaceae [70–72]. The transient accumulation of acetate with highest peak concentration on day 7 followed by its fast methanogenic consumption was closely linked to the increase in both the activity (defined by increased mRNA abundance) of Methanosarcinaceae (days 7 to 14, Fig. ​Fig.4b)4b) and the CH₄ production rate (days 11-21, Fig. ​Fig.1d).1d). Thereafter, the decline in Methanosarcinaceae’s activity (days 14 to 28, Fig. ​Fig.4b)4b) corresponded well to a transient decrease in CH₄ production (days 21 to 28, Fig. ​Fig.11d).

Propionate and butyrate are the next important intermediates of organic matter conversion under anoxic conditions [73]. Acetate, H₂, and CO₂ generated by the syntrophic conversion of propionate and butyrate can feed acetoclastic and hydrogenotrophic methanogens. The oxidation of propionate (G₀' = +76 kJ/mol) and butyrate (G₀' = +48,3 kJ/mol) is highly endergonic under standard conditions but can be accomplished by syntrophy with an H₂-utilizing methanogen, which maintains a low H₂ partial pressure [74]. The rapid consumption of butyrate between days 7 and 11 (Fig. ​(Fig.1b)1b) corresponds well to the initial activity increase of H₂-utilizing Methanocellaceae (Fig. ​(Fig.4b).4b). The net consumption of propionate from day 14 towards day 28 indicates that in addition to butyrate, propionate fueled acetoclastic methanogenesis (Fig. ​(Fig.11c).

Composition and activity dynamics of the methanogenic community

Our study identified three successional phases defined by two activity peaks separated by an intermittent decrease of methanogenic activity (Fig. ​(Fig.4).4). Methanogenic archaea obtain energy for growth by converting C1 and C2 compounds, including CO₂, formate, acetate, ethanol, methanol, and other methylated compounds, to methane [75–77]. Thus, the succession of methanogen guilds is due to changes in available substrates and their utilization for energy conservation. Methanosarcinaceae were the first methanogens to be stimulated in response to substrates (acetate) released by the anaerobic degradation of rice straw and may be considered early- or rapid-responding methanogens (Figs. ​(Figs.1a1a and ​and4b).4b). With the ongoing anaerobic degradation of complex carbon, the activity of Methanocellaceae as intermediate responders was stimulated, primarily due to their role as syntrophic methanogen partners in the bacterial conversion of propionate to acetate, H₂, and CO₂ (Figs. ​(Figs.1c1c and ​and4b)4b) [78–81].

Methanosarcinaceae and Methanotrichaceae are an excellent example of how two family-level groups compete for the same substrate. Methanosarcinaceae prevailed in acetoclastic methanogenesis at the high acetate concentrations during the early phase, but after acetate declined to a very low level, the late-responding Methanotrichaceae significantly increased in competitiveness and acetoclastic activity (Figs. ​(Figs.1a1a and ​and4b).4b). The abundance dynamics of total mRNA affiliated to either Methanosarcinaceae or Methanotrichaceae agreed well with changes in the relative expression level of their pathway genes indicative of acetoclastic methanogenesis (Figs. ​(Figs.4b4b and ​and6b,6b, c). With a higher maximum rate of acetate utilization and maximum growth rate (Y * k), a higher half-saturation coefficient (KS), and a higher yield coefficient compared with Methanotrichaceae, elevated acetate concentrations are favorable for growth and activity of Methanosarcinaceae. By contrast, members of the Methanotrichaceae are superior at low acetate concentrations due to the investment of energy to activate acetate, thereby leading to a lower k and KS [82–84]. Hence, our results are consistent with the concept that both substrate availability and substrate concentration (threshold concept) are key factors controlling the structure and function of methanogenic communities.

A transient activity decline around day 28 was observed in the absolute abundance of biomarkers, including both genes and transcripts of bacterial 16S rRNA and methanogenic mcrA (Fig. ​(Fig.2).2). This activity decline was also well evidenced by changes in the metatranscriptomic abundance of methanogen rRNA and mRNA (Fig. ​(Fig.4).4). Previous research has shown that over the first 4 weeks of straw decomposition, labile polymers are much faster hydrolytically released than recalcitrant compounds; with the latter having a longer residence time [70, 85]. The anaerobic decomposition of labile straw components, such as pectin and xylan, primarily occurs by a bacterial food web attached to rice straw [6, 36]. Exhaustion of the labile straw components leads not only to a change in the microbial colonization of rice straw but also to a transient decline in methanogenic precursors. In particular, hydrolytic and fermentative bacteria are increasingly detached from rice straw and released into the soil [6, 36]. In consequence, the transition from degrading the labile straw components to decomposing recalcitrant C compounds presumably leads to an overall decline in community activity around day 28 until critical sizes of functionally active populations and methanogenic precursor pools are re-established, thereby leading to a second activity peak during the late phase.

Characterization of dominant Methanosarcina populations

The intra-family analysis of Methanosarcinaceae identified four distinct populations, with the Group II population prevailing over the complete 120-day incubation period (Fig. ​(Fig.5).5). This implies the same methanogen population with close affiliation to Methanosarcina sp. MSH10X1 dominated Methanosarcinaceae during the early and late activity phases (Figs. ​(Figs.44 and ​and5).5). An ANI value of 95% is considered species demarcation [65]. Thus, our three Methanosarcina MAGs recovered from slurry samples incubated for 21, 28, and 35 days presumably belong to the same species. Moreover, the high ANI value found between the three Methanosarcina MAGs and Methanosarcina sp. MSH10X1 (around 85%) suggests that the MAGs and strain MSH10X1 represent different but closely related species of the Group II population. (Fig. S6a). Phylogenetic analysis of the 16S rRNA gene located on the Methanosarcina MAG-21 and representative full-length 16S rRNA sequences assembled from our metatranscriptomic datasets further substantiate that the three Methanosarcina MAGs are closely affiliated to Methanosarcina sp. MSH10X1. This provides additional evidence for the dominant role of the Methanosarcina Group II population in the slurry incubations (Fig. S7).

Expression of methanogenic pathways by Methanosarcina spp.

Previous paddy soil studies have demonstrated that acetoclastic methanogenesis is the major methane production pathway operated by members of the Methanosarcinaceae [16, 19, 28]. Indeed, we observed a high expression level of key genes (ack, pta, and cdhA-E) related to acetoclastic methanogenesis by Methanosarcina populations. The expression of ack and pta reached their peak abundance in the early phase, while their transcripts were detected only on a low level in the late phase. Obviously, acetoclastic methanogenesis was the major pathway operated by Methanosarcina spp. in the early phase (Fig. ​(Fig.6).6). The competitive mapping results further confirmed the major role of acetoclastic methanogenesis in Methanosarcina populations (including Groups I, II, and IV) during the early phase (Fig. ​(Fig.77b).

However, we also detected the expression of key genes encoding methylotrophic methanogenesis (mta, mtb, mtm, and mtt), but in particular transcripts of mtaBC (accounting for 0.8 to 5.2% of total mRNA affiliated to the KEGG level 3 category “methane metabolism”). The enzyme methyltransferase/methanol corrinoid protein encoded by mtaBC is a specific biomarker for methanol-dependent methanogenesis (Fig. ​(Fig.6b).6b). The transcripts involved in methylotrophic methanogenesis were affiliated with Methanosarcinaceae. While their greatest expression level agreed largely with the two activity peaks of Methanosarcinaceae, the first peak abundance of methylotrophic mRNA showed a time shift relative to the peak transcript abundance of acetoclastic methanogenesis in the early phase (Figs. ​(Figs.7b7b and ​and4b).4b). Thus, it is reasonable to conclude that Methanosarcinaceae produce CH₄ via acetoclastic and methylotrophic methanogenesis in the Philippine paddy soil under anoxic conditions. The competitive mapping approach resulted in the mapping of nearly all methylotrophic mRNA to the Group II representatives, thereby implying that Methanosarcina Group II was the major player in methylotrophic methanogenesis (Fig. ​(Fig.7b).7b). Detection of the genes encoding methylotrophic methanogenesis on the three Methanosarcina MAGs further substantiates this conclusion (Fig. ​(Fig.7a).7a). Actually, Methanosarcina strain MSH10X1 prefers methanol and trimethylamine as substrates over acetate, H₂/CO₂, and other methylated compounds [86]. The utilization of methanol and trimethylamine resulted in the greatest methane yield produced by strain MSH10X1. This pure-culture-based finding agrees well with our result that Methanosarcina Group II is most competitive under conditions favoring methylotrophic methanogenesis. The expression dynamics of the acetoclastic and methylotrophic pathways by the Group II population strongly differed over the 120-day incubation period, thereby providing evidence that they were differentially regulated in response to major changes in substrate availability (Figs. ​(Figs.66 and ​and77b). In addition to methanol released during pectin degradation, certain amounts of methanol released during the decomposition of xylan [87, 88] and lignin [89–91] may have fed methylotrophic methanogenesis. Recalcitrant lignin may have been the primary source of methanol for methylotrophic methanogenesis during the second expression peak around incubation day 60. Here, it is noteworthy that compared to other major rice farming areas (China, Italy), the Philippine rice field soil is particularly rich in humic acids [32]. The anaerobic decomposition of lignin and the release of methanol requires the activity of bacterial laccases-like enzymes that anaerobically depolymerize lignin into smaller polymer units to which the enzyme pool of the anaerobic bacteria can get access to decompose them into monomers [90, 91]. Among the bacterial community, members of the Geobacteraceae are the most promising candidates for being involved in the anaerobic lignin degradation. These bacteria showed not only the expression of putative laccase genes but also their greatest relative expression level during the later phase (days 28 and 35) (Table S21). Thus, a functional interplay may have occurred between members of the Geobacteraceae and the Methanosarcina Group II population, a view that is further supported by the correspondence of their abundance dynamics on both rRNA and mRNA levels (compare Figs. ​Figs.33 and ​and4;4; and see Additional File 3 for further details). Another putative candidate for being involved in the degradation of recalcitrant lignin may be members of the Peptococcaceae, given the reasonable correspondence of their rRNA and mRNA peak abundances with the expression peaks of methylotrophic methanogenesis. The Peptococcaceae have been identified as putative aromatic ring-cleavage bacteria in rice field soil under lignin-degrading methanogenic conditions [92].

We are grateful to Nicole Paczia and Peter Claus (Core Facility for Metabolomics and Small Molecule Mass Spectrometry, MPI for Terrestrial Microbiology) for their technical support. Xin Li received a stipend from the Chinese Scholarship Council and was a member of the International Max Planck Research School Marburg.