Disentangling Syntrophic Electron Transfer Mechanisms in Methanogenesis Through Electrochemical Stimulation, Omics, and Machine Learning

Background: Interspecies hydrogen transfer (IHT) and direct interspecies electron transfer (DIET) are two syntrophy models for methanogenesis. Their relative importance in methanogenic environments is still unclear. Our recent discovery of a novel species Candidatus Geobacter eutrophica with the genetic potential of IHT and DIET may serve as a model species to address this knowledge gap. Results: To experimentally demonstrate its DIET ability, we performed electrochemical enrichment of Ca. G. eutrophica-dominating communities under 0 and 0.4 V vs. Ag/AgCl based on the presumption that DIET and extracellular electron transfer (EET) share similar metabolic pathways. After three batches of enrichment, acetate accumulated in all reactors, while propionate was detected only in the electrochemical reactors. Four dominant fermentative bacteria were identied in the core population, and metatranscriptomics analysis suggested that they were responsible for the degradation of fructose and ethanol to propionate, propanol, acetate, and H 2 . Geobacter OTU650, which was phylogenetically close to Ca. G. eutrophica, was outcompeted in the control but remained abundant and active under electrochemical stimulation. The results thus conrmed Ca. G. eutrophica’s EET ability. The high-quality draft genome (completeness 99.4%, contamination 0.6%) further showed high phylogenomic similarity with Ca. G. eutrophica, and the genes encoding outer membrane cytochromes and enzymes for hydrogen metabolism were actively expressed. Redundancy analysis and a Bayesian network constructed with the core population predicted the importance of Ca. G. eutrophica-related OTU650 to methane production. The Bayesian network modeling approach was also applied to the genes encoding enzymes for alcohol metabolism, hydrogen metabolism, EET, and methanogenesis. Methane production could not be accurately predicted when the genes for IHT were in silico knocked out, inferring its more important role in methanogenesis. Conclusions: Ca. G. eutrophica is electroactive and simultaneously performs IHT and DIET. The results from the metatranscriptomic analysis have provided valuable information for enrichment


Background
Anaerobic digestion is widely used to convert high-strength waste streams to biogas. This renewable energy source is composed of only ~60% methane (CH 4 ) and needs downstream upgrade to increase the heating value [1]. Therefore, it is of practical signi cance to understand the metabolic pathways that drive biogas production and develop engineering strategies to increase CH 4 content [2]. In anaerobic digestion, CH 4 is produced predominantly via interspecies hydrogen transfer (IHT) and/or an acetoclastic pathway [3]. The former has been extensively studied as a classical model of syntrophy [4], in which hydrogenproducing bacteria use proton as an electron sink for energy conservation, and hydrogenotrophic methanogens scavenge the produced hydrogen gas (H 2 ) for methanogenesis [5]. Of the 155 methanogen cultures, 120 can use H 2 as an electron donor, and 110 are strictly hydrogenotrophic [6], underpinning the importance of IHT in many methanogenic environments.
Direct interspecies electron transfer (DIET) was recently discovered to be a novel syntrophic model of methane production with Geobacter acting as a key syntroph [7,8]. In an early study, the aggregates from an up-ow anaerobic sludge blanket reactor were characterized to be conductive and dominated by Geobacter spp. [7]. Adding conductive media into anaerobic digesters further accelerated methane production [9,10]. Because many Geobacter spp. were capable of extracellular electron transfer (EET) [11][12][13], they were hypothesized to transfer electrons directly to methanogens through direct contact. This mechanism was later proven valid by co-culturing Geobacter metallireducens with Methanosaeta harundinacea or Methanosarcina barkeri [8,14].
The genus Geobacter is found dominant in many digester communities and may play a key role in methane production [15][16][17]. Multiple lines of evidence suggest that Geobacter spp. can form syntrophic associations with methanogens via IHT. For example, Geobacter sulfurreducens can oxidize acetate and produce H 2 in the presence of hydrogen-oxidizing partners [18,19]. This species has been reported to produce H 2 in bioelectrochemical systems [20], leading to the enrichment of hydrogenotrophic methanogen [21]. From the G. sulfurreducens genome, four distinct [NiFe]-hydrogenase-encoding gene clusters that are putatively involved in hydrogen-dependent growth have been identi ed [22].
It is reasonable to hypothesize that the Geobacter spp. dominant in methanogenic environments may transfer electrons simultaneously via IHT and DIET [23,24], which allows them to be ecologically versatile and thrive under varying conditions. This hypothesis can be strengthened by our recent study of anaerobic reactors amended with conductive granular activated carbon (GAC), where we observed high abundances of Geobacter-related taxa (up to 55%) and enhanced methane production [25]. The most dominant population was con rmed to be a novel Geobacter species (Candidatus Geobacter eutrophica strain bin GAC1) based on the phylogenetic and phylogenomic analyses. Its draft genome revealed the genetic potential of both IHT and DIET, including the ech cluster for energy-conserving [NiFe] hydrogenase complex [26], a cytoplasmic bidirectional NAD-reducing HoxPLSUFE complex, a periplasmically oriented membrane-bound HyaSLBP complex [22], and 38 protein-coding sequences for c-type cytochromes that were putatively involved in EET [27].
The extent to which IHT and DIET contribute to methanogenesis remains unclear and is currently evaluated using a predictive model [28,29]. Ca. G. eutrophica may serve as a model species to address this knowledge gap. The objective of this study is to experimentally con rm Ca. G. eutrophica's DIET ability through electrochemical stimulation based on the working hypothesis that it uses similar metabolic pathways for DIET and EET [12,30]. If Ca. G. eutrophica performs DIET, it may also use an anode electrode as an electrode acceptor and can be electrochemically stimulated and enriched. To this end, Ca. G. eutrophica-dominating GAC was collected from our previous study and enriched in twochamber bioelectrochemical systems with the electrode poised at positive potential [25]. After three batches of enrichment, we sequenced the 16S rRNA and rRNA gene to characterize the overall microbial activity and community structure, respectively. We also performed metagenomic and metatranscriptomic analyses to understand the ecophysiology of the dominant Geobacter and their fermentative and methanogenic partners. Finally, we developed a Bayesian network approach to provide predictive insights into the importance of IHT and DIET to methanogenesis.

Reactor setup and operation
Two-chamber bioelectrochemical systems were built as previously described [31]. The anode (total volume 130 mL) was amended with 20 mL of fresh GAC and inoculated with 5 mL of cultivated GAC. The inoculum was collected from an up-ow packed-bed bioreactor that had been operated for three months according to our previous study [25], which resulted in a high abundance of Ca. G. eutrophica. The anode electrode was a stainless-steel brush immersed in GAC. An Ag/AgCl electrode [3M KCl, 0.21 V vs. SHE (standard hydrogen electrode)] was installed adjacent to the anode electrode as a reference. The reactors were fed with 100 mL of substrate mimicking soft drink-processing wastewater [17]. The substrate contained 1500 mg/L fructose and 1100 mg/L polyethylene glycol (PEG) 200 as the carbon sources and was slightly modi ed by adding 50 mM phosphate buffer saline to buffer the pH change. The chemical oxygen demand (COD) of the in uent was 3000 mg/L. The catholyte was 100 mM phosphate buffer saline. The reactors were operated under a batch mode with a hydraulic retention time of three days.
To prepare seeds for enrichment, the anode potential was poised at 0 V vs. Ag/AgCl using a multi-channel potentiostat (MPG-2, BioLogic). The control reactors were operated under an open-circuit condition. Enrichment started after 60 days of incubation (20 cycles), and 5 mL GAC was collected to seed Batch 1. The electrochemically stimulated reactors were further operated under three conditions: 0 V, 0.4 V, and 0 V + ethanol. The COD of the ethanol-based substrate was adjusted to 3000 mg/L. The enrichment process was repeated three times, and each lasted 20 cycles (60 days) to ensure the establishment of a stable microbial community. All reactors were operated in triplicate at 37 ºC. A schematic of the enrichment process is shown in Supplementary Figure S1.

Chemical analysis
Samples were collected in the last three cycles of each batch. Biogas produced by the reactors was collected using water displacement for volume and rate measurement. After moisture removal, CH 4 concentration was determined using a GC-2014 Gas Chromatograph (Shimadzu) equipped with a Molecular Sieve 13X packed column (2000 × 2 mm, Restek) and a thermal conductivity detector. For chemical measurements, e uent samples were ltered through 0.22 µm membranes. Soluble COD was measured with a COD digestion kit and DR/4000 U Spectrophotometer (HACH). Total organic carbon (TOC) was measured with a TOC-V Analyzer (Shimadzu). Volatile fatty acids (VFAs) were measured with a 1200 Series HPLC System equipped with a Hi-Plex H column for organic acids (Agilent). Ethanol concentration was measured with a GC-2014 Gas Chromatograph (Shimadzu) equipped with a ame ionization detector. The pH of the e uent was measured using a benchtop pH meter (Fisher). Current production resulted from poised potential was recorded with the manufacturer's software (BioLogic), and Coulombic E ciency (CE) was calculated as described before [32].

Nucleic acid extraction and sequencing
Triplicate GAC samples were collected at the end of each batch and stored immediately at -80 °C before nucleic acid extraction. DNA extraction, library preparation, amplicon sequencing, and microbial community analyses were performed as previously described [33,34]. RNA was extracted following the manufacture's instruction (Qiagen) and puri ed with both RNase-free DNase Set (Qiagen) and Turbo DNAfree (Ambion). Genomic DNA in the puri ed RNA was not detected using DNA-targeting PCR, and the treated products were reverse transcribed to cDNA using a random hexamer (Invitrogen). The V4-V5 region of the 16S rRNA gene was ampli ed with the prokaryote universal primer pair 515F/909R. Puri ed PCR products were pooled in equalmole ratios and sequenced on an Illumina Miseq Bulk 2×300 bp paired-end platform. One reactor under each operating condition was collected from Batches 1 and 3 for metagenomic sequencing (8 samples), and all three reactors were collected from Batches 1 and 3 for metatranscriptomic sequencing (24 samples). Sequencing was performed on an Illumina NovaSeq 6000 platform to generate 2×250 bp paired-reads. Sequencing data are available in GenBank under the accession number PRJNA689312.

Community analysis
Paired-end sequences were assembled and denoised using QIIME 2, and operational taxonomic units (OTUs) were picked using DADA2 [35,36]. Taxonomy was assigned using the QIIME 2 plugin featureclassi er and the Greengenes database as the reference, with the classi er trained on 97% OTUs [37,38].
Statistical analyses, including two-sample t-test, linear regression, principal coordinate analysis (PCoA, based on weighted UniFrac distance) [39], and redundancy analysis (RDA), were perform using R. PCoA and RDA were performed on relative abundance (16S rRNA gene) and microbial activity (16S rRNA), respectively. Permutational multivariate analysis of variance (PERMANOVA) was carried out to test for the signi cant difference of PCoA results (N = 999). The relationships between environmental parameters and taxonomic data in RDA were tested for signi cance with an Analysis of variance (ANOVA) like permutation test (N = 999). A p-value of <0.05 was used to identify a signi cant difference. A total of 705 OTUs were obtained. Core populations were selected based on the following criteria: average relative abundance >0.5% and occurrence >50% across both DNA and RNA samples. A phylogenetic tree was constructed using the neighbor-joining method provided in the ARB program [40].

Metagenomics and metatranscriptomics
Metagenomic reads were trimmed using Trimmomatic with a quality score of 30, sliding window 6 bp, and a minimum length of 200 bp [41]. To retrieve genomes with high quality, the trimmed reads were processed using a method adapted from a recent study [42]. Brie y, each metagenomic datum was individually assembled, and all data were co-assembled using SPAdes [43]. The assembled contigs were binned using MaxBin2, MetaBAT2, and CONCOCT with default parameters [44][45][46]. The genomes from the three binning strategies were consolidated using metaWRAP, and the consolidated genomes from the two assembly strategies were combined using dRep [47,48], with a quality cutoff of >98% completeness and <2% contamination using CheckM [49]. A nal re-assembly step was performed using SPAdes to further improve the genome quality. Taxonomy classi cation of the metagenome-assemble genomes (MAGs) was performed using CAT [50], and the Geobacter populations were extracted to build a phylogenomic tree by including the draft genome of Ca. G. eutrophica and publicly available Geobacter genomes using Phylophlan [51]. Genes were predicted using Prodigal and annotated using Prokka [52,53]. A Hidden Markov Model database composed of 36 sequences for the conductive type IV pilA gene [54], 7 for the genes encoding triheme periplasmic c-type cytochromes (Ppc), and all the genes for the outer membrane c-type cytochrome family (Omc) were manually curated. The Geobacter genomes were searched against the database using HMMER for identi cation of putative genes involved in EET [55,56]. For gene expression analysis, metatranscriptomic reads were trimmed using Trimmomatic with a quality score 30, sliding window 10 bp, and minimum length 100 bp. Ribosomal RNA was removed using SortMeRNA [57], and the ltered reads were mapped to MAGs using the Bowtie 2 [58]. The gene expression level was calculated for individual bins as reads per kilobase transcript per million reads (RPKM) using featureCounts [59].

Bayesian network analysis
Bayesian networks were constructed as previously described [33]. For the network constructed with 16S rRNA data, abundances of the core population and environmental parameters were combined as a single matrix and normalized to between 0 and 1 using Equation 1: where v_norm ij is the normalized variable j in sample i, v ij is the observed variable j in sample i, min(v j ) and max(v j ) are the minimum and maximum values of variable j. Bayesian networks were built with R package "bnlearn" using hill-climbing algorithm, and the parameters of the networks were calculated with Maximum Likelihood parameter estimation method. The networks were examined with leave-one-out cross-validation [60]. Each cross-validation generated a simpli ed microbial community composed of the core population. Bray-Curtis similarity between the predicted and actual community was calculated and compared with those obtained from a null model. The null model was trained by using average taxa abundance [61]. The correlation between the observed and predicted environmental parameters was performed. The prediction was also validated using relative root-mean square error (RMSE), with y as the experimental values, y max as the maximum experimental value, and t as the predicted values: After validation, Bayesian networks at the OTU and genus levels were constructed from all datasets. To construct networks with metatranscriptomic data, genes encoding proteins for alcohol metabolism, hydrogen metabolism, EET (putative outer membrane c-type cytochromes and conductive pili), and methanogenesis were selected as the input from dominant fermentative bacteria, Geobacter spp. and Methanobacterium spp. The model was built following the knowledge about interspecies interaction obtained from community analysis and omics, and the nodes were directed strictly from fermentative bacteria/Geobacter over methanogens to CH 4 using a blacklist function in the bnlearn package. To statistically infer the importance of DIET and IHT, the associated genes were manually silenced (i.e., the expression level was set to 0).

Results
Effects of enrichment on system performance Before enrichment, a cultivation step was performed to alleviate potential electrochemical shock to the methanogenic population by transferring Ca. G. eutrophica-containing GAC from a packed-bed bioreactor to bioelectrochemical systems amended with fresh GAC. After 20 cycles (60 days) of cultivation, the electrochemical seed (0 V vs. Ag/AgCl) showed lower methane production than the control (0.1 L/L/d vs. 0.5 L/L/d, Figure 1A) but similar COD removal of about 1200 mg/L ( Figure 1B). Additionally, acetate accumulation was observed in the electrochemical seed e uent (up to 50% of the TOC, Figure 1C). The results indicate strong electrochemical shock to the methanogenic population, particularly acetoclastic methanogens, during the initial cultivation. GAC was collected from the seeds for enrichment.
In Batch 1, the methane production rate of the control reactors increased slightly to 0.6 L/L/d ( Figure 1A), while that of the 0-V reactors increased signi cantly to 0.4 L/L/d, suggesting recovery of the methanogenic population from electrochemical shock. One of the 0-V reactors produced low biogas throughout the rst batch due to rapid pH drop (average 5.2 after each cycle, Supplementary Figure S2), but the reason was not clear. The 0-V reactors exhibited more stable performance in the following batches and produced CH 4 at a rate of 0.6 L/L/d. Higher potential (0.4 V vs. Ag/AgCl) did not affect methane production, and both 0-V and 0.4-V reactors produced CH 4 slightly slower than the control (p <0.05). Another group of reactors were operated at 0 V but fed with ethanol as the sole carbon source, as described in pure-culture DIET studies [8]. However, the EtOH reactors produced CH 4 at a much lower rate of 0.1 L/L/d in all batches.
At the early stage of the enrichment, the fructose/PEG-fed reactors (i.e., control, 0 V, and 0.4 V) showed similar COD removal of about 1200 mg/L, whereas the EtOH reactors removed only 250 mg/L COD. COD removal doubled for all reactors in the following batches, but this enhancement was not re ected by CH 4 production, which increased only by <20% ( Figure 1A). In the meantime, CE showed a decreasing trend and dropped by ~3% for the fructose/PEG-fed electrochemical reactors and by 10% for the EtOH reactors. High COD removal and low CE are expected as a signi cant proportion of the electrons will be channeled to biomass, metabolites, and CH 4 when fermentable substrates are fed as the electron donor [62].
The concentrations of short-chain VFAs (acetate, propionate, and butyrate) and ethanol in the e uent were converted to TOC for comparison ( Figure 1C). For the electrochemical reactors, e uent TOC dropped from 800 mg/L in the seeds to about 500 mg/L in Batch 1. The abnormally high TOC found in one of the 0-V reactors in this batch was consistent with its gas production and COD removal and was attributed to the low pH that caused reactor failure. Total VFAs in those reactors further decreased to below 200 mg-TOC/L in Batch 2 and 3 with noticeable propionate accumulation. Ethanol was detected with a low concentration (<10 mg-TOC/L), which was not in agreement with our previous study [25]. For the EtOH reactors, e uent TOC remained stable throughout the enrichment and consisted of 400 mg-TOC/L acetate followed by 100 mg-TOC/L ethanol. Acetate accumulation and low CH 4 production together suggest inhibition of acetoclastic methanogenesis, likely because the presence of acetate and electrode favors the growth of electroactive acetate scavengers and potentially leads to the washout of acetoclastic methanogens [63].

Effects of enrichment on microbial population and overall activity
A clear shift in the community structure in the 0-V and 0.4-V reactors can be seen from the 16S rRNA gene-based PCoA results (Supplementary Figure S3). The communities in those reactors became signi cantly different from the seed communities (PERMANOVA, p <0.05) and clustered closely with the control communities in Batch 3, demonstrating successful enrichment. Similar results were reported by a previous study, in which poised potentials exerted minor impacts on Geobacter-dominating microbial communities [64]. The EtOH-fed communities, on the other hand, were separated from the fructose/PEGfed communities since Batch 1 and remained stable throughout the enrichment. The signi cant difference in microbial communities agreed well with the system performance (Figure 1), highlighting the key role of the substrate as a deterministic factor that drives microbial community dynamics [34,65].
A core population composed of 22 OTUs was selected based on the criteria of average abundance >0.5% and occurrence >50% across all DNA and RNA samples ( Figure 2). In the 0-V and 0.4-V reactors, Ethanoligenens OTU682 was initially abundant (57% 16S rRNA gene abundance) and active (49% 16S rRNA abundance) and was replaced by Clostridium OTU92 during the enrichment. Rhodocyclaceae OTU196 and Desulfovibrio OTU66 were the dominant taxa in the EtOH reactors. G. sulfurreducens-related OTU268 was stimulated and enriched in all electrochemical reactors, and the correlation with CE (Supplementary Figure S4) revealed its electroactive nature [66]. Also enriched were strict hydrogenotrophic Methanobacterium spp. [67,68].
Geobacter OTU650 is of particular interest due to its high phylogenetical similarity to Ca. G. eutrophica (Figure 2). It was outcompeted in the control but not under electrochemical stimulation. Although its abundance decreased to between 5% -10% in the 0-V and 0.4-V reactors in Batch 3, OTU650 was still among the top three most abundant taxa. Moreover, it contributed to 15% to 20% of the community 16S rRNA and showed high activity throughout the enrichment. OTU650 was absent in the EtOH reactors and might not be a competitive utilizer of ethanol and acetate while respiring on an electrode. Instead, Activity-based RDA predicted the association between this taxon and propionate (Supplementary Figure  S4). OTU650 is also predicted to be involved in CH 4 production. The signi cant difference between the control and 0-V/0.4-V reactors (p <0.05) con rms that Ca. G. eutrophica-related OTU650 is capable of EET.

Fermentative bacteria providing substrate for methanogenesis
Metagenomic analysis yielded 28 high-quality MAGs (>98% completeness and <2% contamination, Supplementary Figure S5), whose bin coverage and percentage of mapped reads in Batches 1 and 3 agreed well with the relative abundance of the selected core population shown in Figure 2. We reconstructed the metabolic pathways for fructose utilization, propionate/propanol accumulation, ethanol production/utilization, and H 2 metabolism for the four primary fermentative bacteria (Supplementary Figure S6).
Ethanoligenens co_bin 13 was abundant in the 0-V and 0.4-V reactors in Batch 1 (Supplementary Figure  S5). The genes for fructokinase (scrK), mannose PTS system (manXZ), and mannose isomerase (manA) were highly expressed, suggesting that co_bin 13 could metabolize fructose directly or with mannose as an intermediate. The resulting β-D-fructose 6-phosphate was fed into glycolysis and led to the accumulation of propionate through the generation of glycerone phosphate, methylglyoxal, and propionyl-CoA (Supplementary Figure S6). It actively expressed six genes that are important for the conversion of propionate/propionyl-CoA, but the genes for downstream propionyl-CoA utilization were not detected. Similar to other Ethanoligenens species [69][70][71], co_bin 13 could carry out complete glycolysis to produce ethanol, as indicated by the expression of the adhE gene encoding aldehyde-alcohol dehydrogenase. However, genes for H 2 -evolving hydrogenase were not found in its draft genome. The results collectively suggest that Ethanoligenens co_bin 13 ferments fructose to ethanol, propionate, and propanol.
Clostridium B1_C1_bin17 was ubiquitously abundant in the fructose-fed reactors (Supplementary Figure  S5). In addition to the two fructose metabolism pathways mentioned above, Clostridium B1_C1_bin17 expressed a third route with mannitol as an intermediate. It could convert the produced β-D-fructose 6phosphate to propionate-CoA either via the glycerone phosphate-methylglyoxal route or through lactate metabolism (Supplementary Figure S6), further underpinning its metabolic exibility. B1_C1_bin17 could also be a major propionate producer. The ability to produce 1,3-propanediol and propanol by Clostridium species was observed in B1_C1_bin17 through the expression of the dhaT gene [72,73], a 1,3-propanediol dehydrogenase-like enzyme that potentially catalyzed NADH-dependent propanal reduction to propanol. Finally, B1_C1_bin17 metabolized ethanol by expressing three alcohol dehydrogenase-encoding genes and the genes for the periplasmic [NiFe] and [NiFeSe] hydrogenases. The expression was one order of magnitude higher in the fructose-fed reactors than in the EtOH reactors, suggesting that it could couple fructose degradation to ethanol metabolism with proton as an electron sink. Overall, Clostridium B1_C1_bin17 contributes propionate, propanol, and H 2 to methanogenesis. Rhodocyclaceae B3_EtOH1_bin9 and Desulfovibrio B3_EtOH1_bin5 preferred ethanol as the substrate (Supplementary Figure S5). Genomic analyses con rmed that the Rhodocyclaceae population was incapable of fructose metabolism, and both were de cient in propionate production (Supplementary Figure S6). The Rhodocyclaceae population showed the ability to metabolize ethanol and carried the genes for the NAD-reducing Hox complex and a periplasmic [NiFeSe] hydrogenase, which indicated its role as a major H 2 -donating partner. The Desulfovibrio population also actively expressed the genes for periplasmic [Fe] and [NiFe] hydrogenase complexes that allowed it to use proton as an electron acceptor and grow syntrophically with H 2 -scavenging partners [74,75]. Alternatively, under electrochemical stimulation, Desulfovibrio B3_EtOH1_bin5 may use its own H 2 as an electron donor to respire on the poised electrode. B3_EtOH1_bin5 showed high activity of the Hnd complex that was recently reported to transfer electrons from H 2 to NADH through avin-based bifurcation [76]. We also observed high activity of the menaquinone reductase complex (Qrc) involved in sulfate respiration [77], which might explain the EET ability found in several Desulfovibrio spp. [78,79].
Geobacter playing a key role in methanogenesis Geobacter B1_C1_bin15 (OTU650) was recovered with an almost complete genome (completeness >99.4% and contamination <0.6%) and con rmed to be phylogenomically nearly identical to Ca. G. eutrophica (Figure 3). Metabolic construction reveals a complete set of genes for propionate metabolism, seven copies of the dhaT gene for 1,3-propanediol dehydrogenase for propanol metabolism, and high activity of several types of alcohol and aldehyde dehydrogenases (Supplementary Figure S7). These results are consistent with the previous nding that G. metallireducens grows syntrophically with DIET partners on propionate, propanol, and ethanol [15].
Through mapping the metatranscriptomic reads to a manually curated database, we identi ed three outer membrane c-type cytochromes that were slightly more active (up to 2 folds, p <0.05) in the electrochemical reactors after enrichment (Figure 3). Among them, OmcF appears to be required for the transcription of the gene encoding OmcC that is directly responsible for Fe(III) reduction and potentially GAC-mediated EET/DIET in this study [80]. Three copies of the conductive pili-encoding pilA gene were also detected but did not show a noticeable difference under different conditions. One of the pilA genes, with 72% identity and 100% coverage relative to that found in G. metallireducens, was actively expressed, implying its involvement in EET/DIET.
We also observed the expression of the genes for putative H 2 -evolving hydrogenases (Figure 3). The Hox complex was globally active across batches and conditions, while the periplasmic [NiFe] and [NiFeSe] complexes became slightly less active after enrichment but still showed 2-fold higher expression under electrochemical stimulation. Further inspection of the B1_C1_bin15 genome revealed that it carries a quinone-reactive Ni/Fe hydrogenase and an Hnd complex in which the genes (hndABCD) are 55-81% similar (>98% coverage) to those found in Desulfovibrio fructosovorans. Both were highly active in the electrochemical reactors and could be involved in H 2 uptake couple with EET to electrode [76].
Geobacter co_bin3 is phylogenomically close to Ca. G. eutrophica (Figure 3) and seems to compete with Geobacter B1_C1_bin15 for the same ecological niche based on their interactive changes in abundance (Figure 2 and Supplementary Figure S5). This is supported by the gene expression pro le, which shows lower activity of the omcF gene, ve copies of the pilA gene, and the genes for H 2 -evolving hydrogenases and the Hnd complex under electrochemical stimulation.
Unlike the other two members in this genus, Geobacter B1_C1_bin22 is ubiquitously abundant in both fructose-and ethanol-fed electrochemical reactors. Its phylogenomic similarity to G. sulfurreducens is re ected by the carriage of a number of cytochromes from the Omc family (Figure 3), including four copies of the well-characterized pili-associated OmcS that is highly upregulated during EET to Fe(III) oxide and electrode [11]. This population also showed activity of both H 2 -evolving and H 2 -uptake hydrogenases, with the latter potentially playing a central role in its growth. Similar to G. sulfurreducens, B1_C1_bin22 is not likely to grow on ethanol [81]. It expressed only one gene for putative alcohol dehydrogenase and lacked acetaldehyde dehydrogenase. Therefore, Geobacter B1_C1_bin22 may compete with hydrogenotrophic methanogens on H 2 as an electron donor for EET, thereby suppressing methane production and producing a high CE in the EtOH reactors ( Figure 1).

Hydrogenotrophic methanogens with DIET potential
Although the bioelectrochemical systems were designed to stimulate DIET-capable Geobacter, we observed high abundances of several Methanobacterium spp. with high phylogenomic similarity to known species (Supplementary Figure S5) and retrieved near-complete genomes (>99% completeness and <1% contamination). These strict hydrogenotrophic methanogens use three routes to recycle coenzyme M and coenzyme B at the end of the methanogenic pathway (Supplementary Figure S8). The ferredoxin:CoB-CoM heterodisul de reductase (HdrA1B1C1), a homolog of the HdrABC complex commonly found in most methanogens, became less active in the control. On the other hand, they possess the heterodisul de reductase [NiFe]-hydrogenase complex for avin-based electron bifurcation coupled with ferredoxin reduction and H 2 oxidation [82], and the genes encoding the complex (hdrABC/mvhADG) were highly active in all treatments. Membrane-bound heterodisul de reductase (hdrD) and F420 dehydrogenase (fpoD), which were proposed to participate in extracellular electron uptake [6], were also present but not active after enrichment, indicating their weak roles in methane production. An interesting nding is that the Methanobacterium spp. carry up to seven copies of the mvhB gene that are actively expressed across different batches and conditions, and the encoded polyferredoxins with high iron content have been speculated to participate in electron transfer [83,84]. It is possible that the Methanobacterium spp. use polyferredoxins to shuttle extracellular electrons to the MvhADG/HdrABC complex to complete DIET, but the speci c electron uptake and transfer mechanisms need further investigations.
A predictive understanding of the methanogenic population Bayesian network analysis was performed to understand the microbial interactions and the potential functions of key taxa in a given microbial ecosystem [33]. To investigate the effects of input data type on network training, we selected 17 and 20 OTUs from the DNA and RNA datasets (abundance >0.5% and occurrence >50% in each dataset), respectively. The modeling method was rst validated by reconstructing the core populations. The Bray-Curtis similarities between the predicted and observed communities (0.62 for DNA and 0.64 for RNA) were signi cantly higher than those from a null model (Supplementary Figure S9). Although the predictive power might be compromised by functionally redundant taxa at a high taxonomic resolution [33], the simulation was more accurate than that yielded by arti cial neural networks (constructed to predict acid mine drainage communities) [61]. The Bayesian network approach was further validated by correlating the predicted and observed system performance (Supplementary Table S1). Satisfactory predictions were achieved with an average R 2 >0.61 and RMSE of 0.11. Among the parameters examined, methane production was predicted with the highest accuracy, and the R 2 with the RNA dataset reached 0.94. Other parameters such as acetate and ethanol in the e uent were also better predicted with RNA than with DNA. Overall, RNA was a more robust indicator than DNA with a signi cantly higher average R 2 (0.69 vs. 0.61) and lower RMSE (p <0.05), demonstrating the strong connection between microbial activity and system performance.
After the Bayesian network modeling approach was validated with Bray-Curtis similarities and system performance predictions, a nal network was constructed with RNA at the OTU level ( Figure 4). The inference direction from propionate to Ca. G. eutrophica-related OTU650 implies its potential role as a propionate utilizer, which is consistent with the results from RDA (Supplementary Figure S3) and metabolic reconstruction (Supplementary Figure S7). All three Geobacter taxa in the core population were associated with methane production. OTU650 was assigned with the highest positive coe cient (0.42) followed by OTU28 (0.36), con rming their contribution to methanogenesis via IHT and/or DIET ( Figure  3). OTU268 showed a negative association (-0.12) and potentially competed with hydrogenotrophic methanogens on H 2 . As revealed by the metatranscriptomic pro ling, this G. sulfurreducens-related taxon is incapable of complete oxidation of ethanol but actively expressed hydrogenases for H 2 uptake ( Figure   3). The inference could also explain the less accurate prediction of methane production at a higher taxonomic level (i.e., at the genus level, Supplementary Table S1). In the genus-level network, the three physiologically distinct Geobacter OTUs were combined under the same genus, and their individual impacts on methanogenesis were neutralized, leading to a decrease in the predictive power. The inference presents a signi cant step toward interpreting black-box machine learning models by providing appropriate inputs (i.e., 16S rRNA) for model training [85].

Discussion
IHT is a syntrophic electron transfer mechanism ubiquitously found in many methanogenic environments [6]. On the other hand, DIET was believed to be the dominant mechanism in bioreactors amended with conductive media [86]. Model simulation suggests that electron transfer via IHT may be an order of magnitude slower than via DIET, and up to 33% of methane production can result from DIET [28,29].
However, the relative contribution of IHT and DIET to methane production remains inconclusive in both natural and engineered systems. Ca. G. eutrophica capable of IHT and DIET can serve as a model species for understanding the effects of environmental and operational parameters on the contribution of those two mechanisms to methanogenesis.
To this end, Ca. G. eutrophica needs to be enriched and isolated, and the present study has provided key information for its enrichment. The absence of Ca. G. eutrophica-related OTU650 in the EtOH reactors ( Figure 2) suggests that it prefers carbon sources other than ethanol when carrying out EET/DIET. Potential carbon sources include propionate, propanol, and PEG (Supplementary Figure S4 and S7). It was seen from the metatranscriptomic analysis that the genes for aldehyde dehydrogenase (aldHT) and oxidoreductase (mop) are signi cantly higher in fructose/PEG-fed than ethanol-fed reactors (p <0.05). Based on the difference, we speculated that those enzymes and long-chain-alcohol dehydrogenase are putatively involved in PEG degradation [87]. Anaerobic degradation of PEG has long been observed in Pelobacter spp. [88], a member from the order of Desulfuromonadales to which Geobacter also belongs. Previous studies have consistently detected propionate and acetaldehyde as the by-and end-products during PEG degradation [89,90], and the produced propionate can be directly metabolized by Ca. G. eutrophica (Supplementary Figure S7). Whether Ca. G. eutrophica can degrade PEG and couple PEG/propionate degradation to DIET warrants further inspection.
It is di cult to experimentally determine the contribution of IHT and DIET to methanogenesis without an enriched culture. Therefore, we applied the Bayesian modeling approach to the metatranscriptomic data to gain predictive insight into the importance of those two mechanisms. Using the MAGs of the four fermentative bacteria, three Geobacter spp., and three methanogens discussed above, we extracted the genes related to alcohol oxidation (both ethanol and propanol), H 2 metabolism (both H 2 evolution and uptake), and the putative genes for DIET (e.g., the omc genes in Geobacter and the mvhB genes in Methanobacterium) to train the model ( Figure 5A). The Bayesian network is composed of two components: upstream gene-gene interactions for predicting the expression level of the relevant genes in the methanogens, and a downstream sub-network that links the methanogen genes to methanogenesis.
Such a network structure allows us to impose the inference direction strictly from primary/secondary fermentative bacteria over methanogens to methane production, thereby reconstructing the metabolic interactions among the key partners. A complete network could yield a highly accurate prediction of methane production with R 2 = 0.96 and RMSE = 0.11 ( Figure 5B). To statistically infer the importance of DIET and IHT, we manually silenced relevant genes. The prediction power without the IHT-related genes was signi cantly compromised ( Figure 5B), as evidenced by a lower R 2 value of 0.64 and a higher RMSE of 0.24. Manually setting the expression level of the DIET-related genes to zero slightly decreased the R 2 value of 0.92. This in-silico knockout strategy thus suggests that IHT plays a more critical role than DIET in methane production in the electrochemical reactors.

Conclusions
In this study, methanogenic populations were enriched in two-chamber bioelectrochemical systems under different potential and substrate conditions. System performance became stable after three batches of enrichment. The microbial communities in the fructose-fed reactors, regardless of electrochemical stimulation, shared a high similarity and were signi cantly different from those in the ethanol-fed reactors. Primary fermentative bacteria metabolize fructose and ethanol and provide propionate, propanol, and H 2 as the substrate for Geobacter and methanogens. Geobacter B1_C1_bin15 was phylogenomically identical to Ca. G. eutrophica and genetically capable of propionate/propanol oxidation, EET/DIET, and IHT. It was outcompeted in the control but remained highly abundant and active under electrochemical stimulation. Hence, its EET ability was con rmed. Bayesian networks inferred the more important role of IHT in methanogenesis.

Declarations
Ethics approval and consent to participate Not applicable Consent for publication Not applicable Availability of data and materials Sequencing data are available in GenBank under the accession number PRJNA689312.