Metagenomic Analysis of Anaerobic Microbial Communities Degrading Short-Chain Fatty Acids as Sole Carbon Sources

Analyzing microbial communities using metagenomes is a powerful approach to understand compositional structures and functional connections in anaerobic digestion (AD) microbiomes. Whereas short-read sequencing approaches based on the Illumina platform result in highly fragmented metagenomes, long-read sequencing leads to more contiguous assemblies. To evaluate the performance of a hybrid approach of these two sequencing approaches we compared the metagenome-assembled genomes (MAGs) resulting from five AD microbiome samples. The samples were taken from reactors fed with short-chain fatty acids at different feeding regimes (continuous and discontinuous) and organic loading rates (OLR). Methanothrix showed a high relative abundance at all feeding regimes but was strongly reduced in abundance at higher OLR, when Methanosarcina took over. The bacterial community composition differed strongly between reactors of different feeding regimes and OLRs. However, the functional potential was similar regardless of feeding regime and OLR. The hybrid sequencing approach using Nanopore long-reads and Illumina MiSeq reads improved assembly statistics, including an increase of the N50 value (on average from 32 to 1740 kbp) and an increased length of the longest contig (on average from 94 to 1898 kbp). The hybrid approach also led to a higher share of high-quality MAGs and generated five potentially circular genomes while none were generated using MiSeq-based contigs only. Finally, 27 hybrid MAGs were reconstructed of which 18 represent potentially new species—15 of them bacterial species. During pathway analysis, selected MAGs revealed similar gene patterns of butyrate degradation and might represent new butyrate-degrading bacteria. The demonstrated advantages of adding long reads to metagenomic analyses make the hybrid approach the preferable option when dealing with complex microbiomes.


Introduction
Biogas production by anaerobic digestion (AD) provides a sustainable option to transform organic waste into a renewable energy carrier. Under anoxic conditions, organic matter is converted to short-chain fatty acids (SCFAs) as intermediates, and finally to methane and carbon dioxide, catalyzed by different microbes in four steps: hydrolysis, acidogenesis, acetogenesis, and methanogenesis. The first three steps are facilitated by bacteria, whereas the last step is performed by archaea. There are multiple types of biogas reactors. The continuous stirred tank reactor (CSTR) is the most common type for agricultural biogas plants in Germany. The organic loading rate (OLR) is a key parameter utilize the hydrogenotrophic pathway, which accounts for the remaining portion of methane production resulting from the reduction of carbon dioxide with hydrogen/formate [30,31].
In this study, Illumina-based shotgun metagenome sequencing was used to characterize and compare compositional and overall functional profiles of reactor microbiomes being subjected to different feeding regimes and OLRs. By restricting the substrate to SCFA, we aimed to select microbiomes that are less diverse than common AD microbiomes and facilitate the functional analysis of acetogenesis and methanogenesis as the last two AD steps. To analyze the genetic profile of individual members of the microbial community (particularly the bacterial community), Illumina-based short-read and Oxford Nanoporebased long-read assemblies were compared with respect to the recovery of high-quality metagenome-assembled genomes (MAGs). Read length is a limitation of the Illumina platform but the high quality is an advantage. Oxford Nanopore reads are long, but the quality is much lower with an error rate of up to 15%, compared to the Illumina error rate of 0.473% (median) [32][33][34][35]. By applying a hybrid assembly approach, the disadvantages of both techniques can be compensated [36][37][38], hence we expected that the completeness of the MAGs increases compared to Illumina-based MAGs. Additionally, the identification and analysis of MAGs can help predict important functional groups [39][40][41].

Samples
The biogas reactor experimental design was described by Bonk et al. [1]. In brief, five reactor experiments were conducted (R-1 to R-5) using CSTRs with a working volume of 8 L (R-1 to R-3) or 6 L (R-4 and R-5) (total volume 15 L), which were operated under mesophilic conditions (37 • C), and with a hydraulic retention time of 8 days. Acetate, propionate and butyrate were fed as sole carbon sources in a synthetic medium in a ratio of 0.45:0.1:0.45, based on chemical oxygen demand (COD). This ratio reflects the product spectrum of acidogenic fermentation, in which typically more acetate and butyrate is formed than propionate (see e.g., [42]). Feeding was varied across two parameters: OLR (high and low) and feeding regime (continuous and discontinuous). Reactors R-1 and R-4 were continuously fed, while reactors R-2, R-3, and R-5 received a single feeding pulse per day ( Table 1). All reactors were inoculated from a laboratory scale reactor that was fed with the same substrates. The original inoculum was derived from a lab-scale reactor digesting distillers grains. Table 1. Experimental setup and process parameters according to [1]. For discontinuously fed reactors, measurements were taken just prior to feeding; mean values and standard deviation are reported for the stable operating phase (29 days for R-1, R-2, R-3, and 44 days for R-4, R-5) with daily measurements of pH and methane production and weekly measurements of total SCFA and microbial biomass (see Figure S1). COD: chemical oxygen demand; VS: volatile solids.

R-1 R-2 R-3 R-4 R-5
OLR (gCOD/L/d) 1 All reactors had reached a stable production phase when sampling occurred (see Supplementary File S1 for details on observed process performance of all reactors). The following parameters were recorded: pH value, methane production rate, total SCFA concentration, and microbial biomass. The experimental setup, process parameters, and measured data in the stable production phase are summarized in Table 1.
For metagenome analysis using the Illumina MiSeq platform, 1.5 mL reactor effluents were centrifuged for 2 min at −7 • C and 15,000× g, immediately upon sampling. The supernatant was removed. The sample pellets for subsequent DNA extraction were stored at −20 • C. The DNA extraction followed the instruction (buffer SL2, no enhancer solution) of the NucleoSpin ® Soil Kit (MACHEREY-NAGEL GmbH & Co. KG, Germany). DNA purity was assessed with a NanoDrop One spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and for DNA quantification, the Qubit ™ dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA) was used.

Metagenome Sequencing
Libraries for whole-genome shotgun sequencing were prepared using the Nextera ® XT DNA Library Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer's protocol. Based on the quantification results by Qubit™ dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, USA), the DNA samples were diluted to 0.2 ng/µL.
For indexing, the indexes i7 and i5 were used. Samples were purified using 30 µL AMPureXP beads (AGENCOURT ® AMPURE ® XP Kit) according to the manufacturer's instructions. The purified library was eluted in 50 µL resuspension buffer. Finally, the library was checked by the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) applying the Agilent High Sensitivity DNA Kit and by Qubit™ dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA).
The library pool was generated using the concentration (based on Qubit™ quantification) and the average size of DNA fragments (in bp) per sample (based on the Agilent 2100 Bioanalyzer quantification). An individual volume per library was pipetted to obtain a pool library of 50 nM, which was finally controlled using the Agilent High Sensitivity DNA Kit and Qubit dsDNA HS Assay Kit. Libraries including 5% PhiX DNA were sequenced on the Illumina MiSeq platform. For the reactors R-1 to R-3, the MiSeq Reagent Kit v2 (Illumina) with 500 cycles (2 × 250 bp reads) was used, and for the samples from reactors R-4 and R-5, the MiSeq Reagent Kit v3 (Illumina) with 600 cycles (2 × 300 bp reads) was used.
Extracted DNA was prepared for sequencing on the MinIon Mk1B platform (Oxford Nanopore Technologies, Oxford, UK) using the SQK-LSK108 Ligation Sequencing Kit (Oxford Nanopore Technologies, Oxford, UK) according to the manufacturer's protocol. Libraries were loaded on rec V R9. 4  First, the Illumina raw reads were quality-filtered by Trimmomatic v0.36 [43] to a minimum read length of 35 bp and minimum average quality of 20. Quality control before and after filtering was performed with FastQC [44] and compared using MultiQC [45]. Between 5 and 25 million high quality reads per sample were generated.
Basecalling of Nanopore raw reads was performed with Albacore v2.3.1 using default settings. Adapters were removed with Porechop v0.2.3 (https://github.com/rrwick/ Porechop) using an identity threshold of 90% and splitting reads when adapters were found with 90% identity in the middle of a read. Quality control was performed with FastQC and MultiQC as well.

Analysis of Community Composition and Functional Potential
Community composition was analyzed using a pipeline based on information of expected species of the AD process, as previously described by Becker et al. [46]. Expected species were chosen based on a literature search in the context of AD, gene prediction results using Metaxa v2.1.3 [47,48], and amplicon sequencing results from Bonk et al. [1], and combined in a reference database. The FASTA files including genomes of these expected species were included in this reference database [46]. According to the database, the reads were classified into expected species reads (aligned against the database) and unexpected species reads (unaligned). Bowtie2 v2.3.3.1 [49] was applied to align the Illumina reads against the database of expected species for the classification. Expected reads were directly annotated using DIAMOND v0.9.8 [50] and MEGAN6 v6.12.0 [51] based on the NCBI-NR database (accessed on 11 July 2018), whereas unexpected species reads were additionally subjected to de novo assembly and binning steps.
The EggNOG database [58] option of DIAMOND/MEGAN5 was used to annotate the functional potential of expected reads and assembled unexpected reads by using the protein accession mapping file from October 2016. Finally, the read counts of annotations were combined to cover the whole functional potential (see Supplementary File S2, B.2.2).
To generate a functional catalog of all reads per sample, interleaved paired and unpaired reads were annotated together, without differentiating between expected and unexpected species reads as before. The databases used for the annotation with DIAMOND and MEGAN5 were the protein accession mapping file of EggNOG [58] from October 2016 and the GI to NCBI taxonomy mapping files for the NCBI-NR protein database from June 2018 (both available under http://ab.inf.uni-tuebingen.de/data/software/megan6 /download/welcome.html, last accessed: 16 July 2018). Additionally, the taxonomy (genus level) and the functional potential (high-level nodes of EggNOG and Clusters of Orthologous Groups of proteins (COG) IDs) were combined to assess the functional potential per genus (R scripts available at https://git.ufz.de/UMBSysBio/mpp).
Based on the study by Sikora et al. [31], selected enzymes (genes) responsible for the relevant pathways of syntrophic butyrate, propionate, and acetate oxidation as well as methanogenesis, based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, were collected (see Supplementary File S1, Tables S1 and S2). This list of genes was used to compile the functional catalog per pathway of interest (Supplementary File S2, B.5).
Finally, the HUMAnN 2.0 pipeline [73] was used to test for the presence of selected pathways [74] in the Prodigal-annotated hybrid MAGs. Coding sequences were mapped against the UniRef50 EC filtered database using DIAMOND. Default parameter settings were used except for lowering the translated query coverage threshold from 90.0 to 50.0 to account for coding sequences being longer than sequencing reads, which are the typical input for HUMAnN 2.0. Identified gene families were mapped to EggNOG COGs using the humann2_regroup_table script.

Impact of Feeding Regime and OLR on Microbial Community Composition
After quality-filtering, between 46% and 67% of all reads could be aligned to the database of expected species. Assembly and binning of unaligned reads led to eleven bins with a contamination below 5% and a completeness higher than 95%. After reassembly, the average total assembly length (mean over the total number of contigs over all samples) was around 26 Mb, the longest contig was around 347 kb, and the highest N50 value was 7.5 kb over all samples (see Supplementary File S2 for details).
The microbial communities in all reactors were dominated by archaea (around 70% relative abundance, see Figure 1a). R-1, R-2, and R-3 revealed similar archaeal community compositions with comparable relative abundances (see Figure 1c). Methanothrix had a relative abundance of around 75% and Methanoculleus of at least 6%, whereas Methanosarcina was less abundant (around 1%). R-4 showed differences in terms of a lower relative abundance of Methanothrix (41%) and a higher relative abundance of Methanoculleus (20%) compared to R-1, R-2, and R-3. Methanosarina increased in abundance up to around 3% in R-4. More conspicuous was the higher relative abundance of Methanosarcina in R-5, with around 32%.
R-4 and R-5 showed higher proportions of bacteria (around 33%) than R-1, R-2, and R-3 (around 26%, Figure 1b). Nine bacterial phyla constituted the bacterial community, with Firmicutes dominating in all reactors (between around 7% and 13%). Bacteroidetes were more abundant in R-1, R-2, and R-3 compared to R-4 and R-5, with around 6.5% compared to 3.8% relative abundance. Besides Spirochaetes, Cloacimonetes were more abundant in R-4 and R-5 than in R-1, R-2, and R-3. Comparing the bacterial communities on the genus level (Figure 1c,d), Desulfovibrio only contributed significantly to the bacterial community composition in discontinuously fed reactors with low OLR (R-2 and R-3). Mesotoga reached up to around 12% but only under continuous feeding at high OLR. At low OLR under continuous feeding (R-1), Syntrophomonas was the dominant genus, while at discontinuous feeding or higher OLR (R-2 to R-5), its relative abundance was lower. Under low OLR, either one or two bacterial genera contributed more than 30% of the bacterial composition (Syntrophomonas under continuous feeding, and either unclassified Bacteroidetes or Desulfovibrio with Endomicrobium under discontinuous feeding), while at high OLR, the distribution was more even. For more information regarding other taxonomic levels, see Supplementary File S2, B.2.1. Relative abundances were calculated using read counts per taxon. In all panels, not assigned refers to reads without hits on the respective taxonomic level as well as unclassified reads. Other refers to reads assigned to other domains (a), phyla, respectively, genera, with relative abundances below 1% (b,c), and bacterial genera with relative abundances below 5% (d). Panels (c,d) share the same legend.

Combining Taxonomic and Functional Annotation
Compared to PCR-dependent analysis methods, metagenomes enable the analysis of community composition without PCR bias, and additionally reveal the functional potential. Even though the bacterial community composition was different in the five reactors, the functional potential was very similar (see Supplementary File S2, B.2.2). Therefore, annotation results were merged over all reactors to assess the functional potential of individual taxa by combining functional and taxonomic annotation of metagenomic reads. For this analysis, 7000 genes and 3500 taxa were considered (see Supplementary File S2, B.4). As SCFAs were supplied as the sole carbon source, the further analysis focused on syntrophic SCFA oxidation and methanogenesis pathways (see Figure  2 and Supplementary File S2, B.5). Relative abundances were calculated using read counts per taxon. In all panels, not assigned refers to reads without hits on the respective taxonomic level as well as unclassified reads. Other refers to reads assigned to other domains (a), phyla, respectively, genera, with relative abundances below 1% (b,c), and bacterial genera with relative abundances below 5% (d). Panels (c,d) share the same legend.

Combining Taxonomic and Functional Annotation
Compared to PCR-dependent analysis methods, metagenomes enable the analysis of community composition without PCR bias, and additionally reveal the functional potential. Even though the bacterial community composition was different in the five reactors, the functional potential was very similar (see Supplementary File S2, B.2.2). Therefore, annotation results were merged over all reactors to assess the functional potential of individual taxa by combining functional and taxonomic annotation of metagenomic reads. For this analysis, 7000 genes and 3500 taxa were considered (see Supplementary File S2, B.4). As SCFAs were supplied as the sole carbon source, the further analysis focused on syntrophic SCFA oxidation and methanogenesis pathways (see Figure 2 and Supplementary File S2, B.5).

Figure 2.
Taxonomic and functional annotation of metagenome reads. Rows contain pathwayspecific genes for propionate, butyrate, and acetate degradation as well as methanogenesis pathways. The color gradient represents read abundance (read counts per genus and potential function as absolute values). Genera with most matches for the selected pathways are shown. IDs in bold indicate alternatives to searched gene IDs, or were not available (see Supplementary File S1).

Figure 2.
Taxonomic and functional annotation of metagenome reads. Rows contain pathwayspecific genes for propionate, butyrate, and acetate degradation as well as methanogenesis pathways. The color gradient represents read abundance (read counts per genus and potential function as absolute values). Genera with most matches for the selected pathways are shown. IDs in bold indicate alternatives to searched gene IDs, or were not available (see Supplementary File S1). Genera indicated in blue were associated with SCFA oxidation pathways and genera indicated in red with methanogenesis pathways. The "unknown" section represents all reads that were not assigned to any known organism in the NCBI database.
Nine genes selected for the oxidative carbon monoxide dehydrogenase/acetyl-CoA synthase (CODH/ACS) pathway (Supplementary File S1) of acetate oxidation were detected and seven of nine genes were unique for the CODH/ACS pathway. Only Acetomicrobium covered the whole CODH/ACS pathway. A high coverage of this pathway was shown for Clostridium and Desulfotomaculum. The acetate oxidation pathway coupled to the reduction of fumarate (Supplementary File S1) was almost completely covered by Burkholderia, Corynebacterium, and Geobacter.
Except for the missing COG2057 gene (acetate CoA-transferase beta subunit), Syntrophomonas and Clostridium covered all selected genes for butyrate oxidation (Supplementary File S1). However, many other genera covered almost the whole butyrate oxidation pathway as well, for example, Bacillus, Desulfotomaculum, Corynebacterium, and Faecalibacterium, among others.
None of the bacterial taxa covered the whole propionate oxidation pathway (Supplementary File S1). The bacterial genera Bacillus, Clostridium, Desulfotomaculum, and Pelotomaculum covered a great part of the propionate oxidation pathway to a lesser extent (around eight of ten selected genes). Accordingly, expected SPOB like Smithella and Syntrophobacter were rare as reflected by their low relative abundance in the whole community.
Syntrophy-specific functional domains of butyrate or propionate oxidation, such as the extra-cytoplasmic formate dehydrogenase (FDH), FDH accessory protein, membraneintegrated proteins, and ribonuclease (Supplementary File S1) were completely encoded in Syntrophobacter, Syntrophomonas, and Syntrophus as well as in Bacillus and Desulfotomaculum.
In our analysis, Clostridium showed a high functional potential and gene coverage in all pathways of acetogenesis ( Figure 2). This is in line with Clostridium species having been reported as versatile contributors in the processes of AD [75,76].
The hydrogenotrophic methanogenesis pathway was completely covered by Methanosarcina, Methanoculleus, Methanothrix, Methanolobus, Methanohalophilus, and Methanospirillum ( Figure 2). The acetotrophic methanogenesis pathway was almost completely covered by Methanosarcina and partially by Methanothrix, Methanohalophilus, and Methanoculleus. The acetate-CoA ligase (COG0365) was encoded in all archaeal genera. Additionally, Methanosarcina also covered the methylotrophic methanogenesis pathway. Therefore, the functional potential of Methanosarcina for all pathways of methanogenesis was detected.

MAGs Reconstructed from Illumina Reads Only
MAGs were reconstructed and analyzed to elucidate the functional role of the detected species in more detail. In total, 94 MAGs were obtained. After refinement, 27 high-quality and 29 medium-quality MAGs were generated (see Table 2), from which 45 bacterial MAGs and 11 archaeal MAGs were allocated (see Supplementary File S2, B.6.1 for details). The archaeal MAGs were mainly assigned to the genera Methanothrix, Methanoculleus, and Methanospirillum. The identified bacterial MAGs were taxonomically more diverse and assigned to several already described MAGs (listed in the GTDB). All results of MAG reconstruction from Illumina MiSeq reads are listed in detail in Supplementary File S2, B.6.1. The mean N50 value was 32.24 kb, the total assembly length around 2.29 Mb, and the largest contig was 93.70 kb.

MAGs Based on the Hybrid Assembly Approach
The Illumina MAG analysis was extended by including Oxford Nanopore reads for the samples from the continuously fed reactors (R-1 and R-4, see Supplementary File S2, B.6.2). In total, 27 MAGs based on the hybrid approach were obtained for both reactors, including 13 high-quality and four medium-quality MAGs. Thirteen of these MAGs could be matched to known species (Table 3). Three high-quality archaeal MAGs were obtained with a completeness between 97.22% and 99.35%, and 0.65% contamination ( Table 3). Two of them represented Methanothrix soehngenii and the third one Methanoculleus. Additionally, ten high-quality bacterial MAGs were reconstructed, see Supplementary File S2, B.6.2. A Cloacimonas acidaminovorans MAG had the highest quality and was circular. Six of the 13 dereplicated MAGs showed a close relationship to Syntrophobacteraceae or Syntrophomonas. Four MAGs represented circular genomes of Methanothrix soehngenii, Syntrophomonas wolfei, Cloacimonas acidaminovorans, and Mesotoga infera B. The mean N50 value was 1.7 Mb, the total assembly length around 3 Mb, and the largest contig was 1.9 Mb.

Pathways of Interest Encoded by MAGs
To elucidate the functional role of organisms represented by the reconstructed MAGs within the community, individual MAGs were tested for the presence of genes encoding syntrophic SCFA degradation and methanogenesis pathways (Figure 3). Four clusters of similar enzyme repertoires were identified. Both members of the first cluster (UFZ-4H3 and UFZ-1H5) were assigned to the family Syntrophobacteraceae. They contained genes of all SCFA oxidation pathways including the almost complete pathways for syntrophic butyrate and propionate oxidation. Two MAGs (UFZ-1H1 and UFZ-4H12), both identified as Methanothrix soehngenii, formed the second cluster and covered a great part of the hydrogenotrophic methanogenesis pathway (eight out of nine genes). Both MAGs contained six of the seven unique genes for hydrogenotrophic methanogenesis (encoding formylmethanofuran dehydrogenase, formylmethanofuran-H 4 MPT formyltransferase, methenyl-H 4 MPT cyclohydrolase, H 2 -forming methylene-H 4 MPT dehydrogenase, and F 420 -dependent methylene-H 4 MPT reductase). Hydrogenotrophic methanogenesis was fully covered in the third cluster (UFZ-1H3 and UFZ-4H8) with both MAGs being assigned to Methanoculleus. The last cluster contained five MAGs (UFZ-4H4, UFZ-1H9, UFZ-1H11, UFZ-4H11, and UFZ-1H7) that were all assigned to the family Syntrophomonadaceae. All enzymes belonging to butyrate degradation were found to be encoded in these MAGs (Figure 3).

Discussion
By using metagenome approaches, we analyzed composition and functional potential of microbial communities in anaerobic digesters under two different feeding regimes at high and low OLRs, with a focus on the conversion of SCFAs to methane. We focused the analysis on selected pathways of particular interest (SCFA oxidation and methanogenesis) and detected associated genes in reads and MAGs. MAGs were reconstructed and annotated to shed light on uncharacterized taxa. The goal was to characterize the full community and identify the putative functional role of not yet described species.
Previous analysis of the methanogenic reactor communities by mcrA-gene based terminal restriction fragment length polymorphism (T-RFLP) analysis had identified Methanothrix and Methanomicrobiaceae as the main constituents, and an increased share of Methanosarcina under discontinuous feeding [1]. These findings were confirmed with metagenome-based analysis in which Methanothrix and Methanoculleus were the main archaeal genera, and the share of Methanosarcina was increased under discontinuous feeding (R-2, R-3, R-5, Figure 1c), a shift that was more pronounced under high OLRs (R-5, [1]), where Methansarcina partially replaced Methanothrix. It is known that Methanothrix is more sensitive to changing environmental conditions and stress conditions (feeding regime or OLR) than Methanosarcina and Methanoculleus [5,[77][78][79][80]. Hence, stress tolerance of AD systems can be promoted if Methanosarcina can be enriched in the reactor at the expense of the more sensitive species, increasing the functional resilience of the methane production process [1]. Discontinuous feeding was shown to increase the relative abundance of Methanosarcina, and Bonk et al. [1] hypothesized that the dynamic environment established by discontinuous feeding provides temporal niches that enable the observed shift in the archaeal community.
The bacterial communities were unique for all reactors, even for reactors R-2 and R-3, which were biological replicates, and no clear trends could be observed regarding OLR and feeding regime (Figure 1d). This has been reported before, including similar values for diversity and evenness of the bacterial communities based on 16S rRNA gene amplicon sequencing data [1].
Using metagenomic data to assess community composition circumvents any PCR biases and, by targeting bacteria and archaea simultaneously, enables the assessment of the ratio between both domains. The observed dominance of archaea was expected as the direct supply of SCFAs as substrate led to the establishment of only the last two steps of AD (i.e., syntrophic SCFA oxidation and methanogenesis), excluding the hydrolysis and acidogenesis steps driven by bacteria. Comparable ratios of bacteria to archaea were also found in other anaerobic digesters supplied with SCFAs [81]. Propionate and butyrate are degraded by bacteria in syntrophy with hydrogenotrophic methanogens, while acetate can either be directly utilized by acetoclastic methanogens or by acetate-oxidizing bacteria, again requiring a hydrogenotrophic methanogen as a syntrophic partner.
Despite the observed variability in bacterial community composition, the overall reactor productivity was not affected, indicating the presence of functional redundancy in the bacterial communities. In low OLR reactors, the bacterial community was dominated by one or two bacterial genera, while high OLR reactors showed a more even distribution of the different genera, indicating that competition played a bigger role under nutrientlimiting conditions. An adaptation of the bacterial community due to an increase in OLR was also observed by Maus et al. [82].
Although being dominated by archaeal genera, the reactors in this study harboured diverse bacterial genera, many with relative abundances below 1%. Given that SCFAs were provided as the sole carbon sources, one could expect that most of these genera are involved in syntrophic SCFA oxidation. Syntrophobacter and Pelotomaculum were identified as known syntrophic oxidizers for propionate. Synthrophomonas was the only butyrateoxidizing genus detected. Genera such as Cloacimonetes [83][84][85], Cryptanaerobacter [86], and Desulfovibrio [87] have been associated with the oxidation of propionic or butyric acid. Candidatus Cloacamonas acidaminovorans, for example, contains all the genes for propionic acid oxidation via methylmalonyl-CoA, but has not yet been isolated [83].
Besides these known and putative SCFA-oxidizers, our MAGs analysis uncovered many bacteria that were apparently not involved in SCFA oxidation. Similar observations have been made in AD reactors being fed with a single SCFA (acetate or propionate) [88].
Yeast extract was present in the medium in the latter study, potentially supporting the growth of other bacteria than SCFA oxidizers. However, this supplement was absent in our study, and hence cannot explain the observed bacterial diversity. Bacteria not directly utilizing the primary carbon source could either be involved in biomass turnover as scavengers or thrive on proteins, amino acids, polysaccharides, and lipids produced by primary consumers and associated syntrophic partners. Autotrophic methanogens can actively excrete amino acids into the medium [75,76], and amino acid exchange has been postulated in a synthetic coculture [89]. Furthermore, metabolite cross-feeding, especially between phylogenetically distant species, might be an important factor driving microbial community dynamics [90], and explain the observed diversity of bacterial taxa, including proteolytic fermenters and amino acid oxidizers not involved in SCFA oxidation. The primary consumers together with their syntrophic partners may hence be the source of excreted nutrients, which sustain a low abundance shadow microbial community of taxa not involved in the primary process.
Overall, the microbial community composition at the genus level was more similar for reactors only differing in feeding regime, whereas those with different OLR harboured more distinct communities, indicating the importance of substrate quantity for shaping community composition. For example, the ratio of bacteria to archaea was mostly determined by the OLR but not by the feeding regime.
The functional potential was similar between reactors despite different community composition. Hence, functional redundancy in SCFA degradation by known or unknown bacterial species could play an important role. By applying a higher OLR (reactors R-4 and R-5), the relative abundance of bacteria increased compared to the low OLR reactors R-1 to R-3. However, the genus Syntrophomonas, comprising typical butyrate-oxidizing bacteria [20,91], showed the highest abundance in the low OLR reactor R-1. Syntrophomonas was also described as the major syntrophic butyrate-oxidizing bacterium in a butyrate-fed system [92]. Its lower abundance in the other reactors indicates that other species took over its function there. Other species only appeared under high OLR conditions or increased their abundance. For instance, Corynebacterium was not found in R-2 and R-3 but had an abundance in the bacterial community of around 13% in R-5, and Mesotoga was present with a relative abundance of less than 1% in R-2 and R-3 but reached around 3% in R-5. Genera like unclassified Bacteroidetes, Candidatus Cloacimonas, Mesotoga, Desulfovibrio, or Endomicrobium were detected with higher abundance, but none of these taxa covered the complete butyrate degradation pathway. Pathway analysis showed that Desulfosporosinus, Clostridium, or Bacillus showed the same pathway coverage of selected genes as Syntrophomonas. However, the relative abundance of these genera was below 1% of the whole community. Thus, potential functional redundancy between yet undescribed taxa and known taxa might be present in our reactors. Candidatus Cloacimonas was detected mainly in the R-4 community with around 5%. Being frequently detected in AD processes [83,93], the reconstruction of a high quality MAG classified as Candidatus Cloacimonas and with 100% completeness (UFZ-4H1) provides the opportunity to identify its role in the process. However, regarding the main pathways analyzed in this study, there was no clear indication for its participation in syntrophic SCFA oxidation. For example, only three genes out of ten for propionate oxidation via the methonylmalonyl-CoA pathway were found (Figure 3), although the complete pathway was reported for the genome of Candidatus C. aminoacidivorans [83].
Regarding typical propionate-degrading bacteria like Desulfotomaculum, Pelotomaculum, Smithella, and Syntrophobacter [20], only Syntrophobacter was detected with more than 5% in the bacterial community and only around 1.5% in the whole community of R-4 and R-5. The pathway analysis of selected genes of propionate oxidation showed a low coverage of these genes by Syntrophobacter. In contrast, Clostridium and Desulfotomaculum displayed a higher share of the typical propionate oxidation genes. A closer examination of Desulfotomaculum suggested that this genus is metabolically versatile, because it encodes many selected genes for acetate, butyrate and propionate oxidation. A single species of Desulfotomaculum has been described as syntrophic propionate oxidizer yet [94].
The results of Illumina reads as well as the analysis of MAGs showed that Methanothrix (UFZ-4H12) encoded all selected genes of the hydrogenotrophic methanogenesis pathway, besides the acetoclastic methanogenesis pathway. This finding was described already by Smith & Ingram-Smith [95]. As hydrogenotrophic methanogenesis activity of Methanothrix has not been described, it has been speculated that this genus can produce methane based on direct interspecies electron transfer [96]. However, as acetic acid was available in our reactors, Methanothrix most likely utilized the acetoclastic methanogenesis pathway.
Read length is a limitation of the Illumina platform but the high quality is an advantage. Oxford Nanopore reads are long but the quality with an error rate up to 15% is much higher compared to the Illumina read quality [32][33][34]. By applying a hybrid assembly approach, the disadvantages of both techniques can be compensated. Even though Nanopore reads get more accurate due to improvements in technology, Illumina, for polishing of especially low-coverage MAGs, is advantageous [97].
We were expecting that the completeness of the genomes would increase compared to Illumina short read-based MAGs. Furthermore, as expected and as described in the literature e.g., [38,97] the contiguity of the MAGs was considerably increased by incorporating long reads.
By the reconstruction of MAGs, several potentially new species could be identified. Focusing on high quality MAGs, a Syntrophobacteraceae species potentially involved in syntrophic SCFA oxidation (UFZ-4H3), three Syntrophomonadaceae species potentially involved in butyrate oxidation (UFZ-4H4, UFZ-1H11, and UFZ-1H7), and one Methanoculleus species potentially involved in hydrogenotrophic methanogenesis (UFZ-1H3) were identified ( Figure 3). The metabolic pathway analysis focused on the two key processes in the reactor: syntrophic SCFA oxidation as the utilization step of the provided carbon sources and methanogenesis as the syntrophically coupled process, providing the required electron sink. However, some of the high-quality MAGs could not be clearly assigned to these key processes (UFZ-4H14, UFZ-1H10, UFZ-4H10, UFZ-4H1, and UFZ-4H6; Figure 3). This indicates that the core community catalyzing the key processes is able to support a wider food web, explaining the observed high diversity. The discovered MAGs provide a first glimpse into this food web, potentially comprising metabolite cross-feeding and scavengers involved in biomass turnover.
The high relative abundance of the archaeal populations limited the read numbers for bacterial populations, limiting the opportunity to reconstruct complete bacterial genomes of rare species. A higher number of bacterial reads would allow to reconstruct more complete genomes of potential new functional groups or species. Nevertheless, 13 medium-and high-quality MAGs could be obtained. The GTDBbk database classification leads to the assumption that some bacterial MAGs belong to the families Syntrophobacteraceae (UFZ-4H3, UFZ-4H4) and Cloacimonadaceae (UFZ-4H1) or the order Bacteroidales (UFZ-1H10).
Concluding, the hybrid assembly revealed to be the optimal approach for genome reconstruction, because the advantage of quality (Illumina approach) and read length (Nanopore) and disadvantages of short reads (Illumina) and lower quality (Nanopore) complement each other. Illumina-based sequencing approaches usually lead to highly fragmented MAGs. This means that complex genomes can only rarely be reconstructed. The functional understanding of microbiomes from AD systems could be expanded and more information could be retrieved from metagenomic samples. After all, an approach that uses a dereplicated set of recovered MAGs from hybrid assemblies is likely to improve our ability to reconstruct genomes and can even enable the strain-specific reconstruction of genomes from metagenomic data.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/microorganisms11020420/s1: Supplementary File S1 (PDF file) containing details on experimental design, reactor performance and enzyme-pathway mappings and Supplementary File S2 (Excel file) containing additional sequencing analysis data.
Author Contributions: D.B. performed DNA extraction, Illumina sequencing and data analysis, and drafted the manuscript. D.P. performed Nanopore sequencing and contributed to data analysis, interpretation of results and revision of the manuscript. F.B. provided the experimental setup and reactor samples and contributed to the revision of the manuscript. S.K. and H.H. contributed to the interpretation of the results and revision of the manuscript. F.C. designed and supervised the study, performed the pathway analysis and contributed to revision of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the German Federal Ministry of Education and Research, grant number FKZ 031A317 (e:Bio project "McBiogas"). Cloud computing facilities used for the bioinformatic analyses were provided by the BMBF-funded de.NBI Cloud within the German Network for Bioinformatics Infrastructure (de.NBI) (031A537B, 031A533A, 031A538A, 031A533B, 031A535A, 031A537C, 031A534A, 031A532B).

Data Availability Statement:
The sequence data have been deposited in the ENA database under BioProject PRJEB42791 (https://www.ebi.ac.uk/ena/data/view/PRJEB42791).