Genome-Centered Metagenomics Analysis Reveals the Microbial Interactions of a Syntrophic Consortium during Methane Generation in a Decentralized Wastewater Treatment System

The application of anaerobic digestors to decentralized wastewater treatment systems (DWTS) has gained momentum worldwide due to their ease of operation, high efficiency, and ability to recycle wastewater. However, the microbial mechanisms responsible for the high efficiency and ability of DWTS to recycle wastewater are still unclear. In this study, the microbial community structure and function of two different anaerobic bioreactors (a primary sludge digestor, PSD, and anaerobic membrane bioreactor, AnMBR) of a DWTS located in Germany was investigated using 16S rRNA gene amplicon and metagenomic sequencing, respectively. The results showed that the microbial community structure was remarkably different in PSD and AnMBR. Methanobacteriaceae and Syntrophaceae were identified as the families that significantly differed in abundance between these two bioreactors. We also used genome-centered metagenomics to predict the microbial interactions and methane-generating pathway, which yielded 21 near-complete assembled genomes (MAGs) (average completeness of 93.0% and contamination of 2.9%). These MAGs together represented the majority of the microbial community. MAGs affiliated with methanogenic archaea, including Methanobacterium sp., Methanomicrobiales archaea, Methanomassiliicoccales archaea, and Methanosaeta concilii, were recruited, along with other syntrophic bacterial MAGs associated with anaerobic digestion. Key genes encoding enzymes involved in specific carbohydrate-active and methanogenic pathways in MAGs were identified to illustrate the microbial functions and interactions that occur during anaerobic digestion in the wastewater treatment. From the MAG information, it was predicted that bacteria affiliated with Bacteroidetes, Prolixibacteraceae, and Synergistaceae were the key bacteria involved in anaerobic digestion. In the methane production step, Methanobacterium sp. performed hydrogenotrophic methanogenesis, which reduced carbon dioxide to methane with hydrogen as the primary electron donor. Taken together, our findings provide a clear understanding of the methane-generating pathways and highlight the syntrophic interactions that occur during anaerobic digestion in DWTS.


Description of a Decentralized Wastewater Treatment System
The DWTS analyzed in this study was built in 2009 in Germany. The main objective of this project was to optimize the processes used for decentralized urban wastewater management. The system was designed to treat wastewater from approximately 175 inhabitants. The detailed process flow of the DWTS is presented in Figure 1. During wastewater treatment, sewage was collected from the vacuum station and pumped into an equalization tank, followed by a sedimentation tank, where the solids were settled for separating treatment. The solids, which were up to 1% to 2% of the influent quantity, were treated in a high-load PSD at approximately 37 °C. To make the treatment more efficient, water was removed during the digestion process by a rotating disk filter. Then, water without settleable solids flowed from the sedimentation tank into an AnMBR. Next, the mixed sludge circulated through an external rotating disk filter, which continuously withdrew the effluent from the bioreactor. Fouling control was achieved by rotating the shaft on which the filtration disks were fixed to create a centrifugal force, which caused the solids to flow off. Finally, the treatment divided the wastewater into a solids-free filtrate, biogas from the two anaerobic bioreactors, and stabilized solids after anaerobic digestion ( Figure 1).

Sampling and Physicochemical Analysis
Water samples were taken from the influent and effluent of the PSD and AnMBR in the DWTS. The physicochemical characteristics were measured as described previously [22]. In brief, influent and effluent chemical analyses (COD, nitrogen species, phosphorus species and alkalinity) were conducted in accordance with the Standard Methods for the Examination of Water and Wastewater, APHA [23]. The pH of the reactors was determined using a DELTA 320 (Mettler Toledo, Columbus, OH, USA). Gas samples were taken from the reactors via a small sampling port located at the top of the reactor. The biogas production was measured with a wet gas meter (COMBIMASS @ GA-s hybrid premium, Binder Group, Germany, and the CH4 production was determined using a gas chromatograph (Clarus 580 Arnel, PerkinElmer, Waltham, MA, USA) equipped with a thermal conductivity detector. Standard methods were used to measure TS, VS, TSS, and VSS [23]. Sludge samples were collected from the AnMBR and PSD and centrifuged at 10000× g for 10 min at 4 °C. The remaining pellet of the sludge was transported on ice to the lab, where it was stored at −70 °C before DNA extraction.

Sampling and Physicochemical Analysis
Water samples were taken from the influent and effluent of the PSD and AnMBR in the DWTS. The physicochemical characteristics were measured as described previously [22]. In brief, influent and effluent chemical analyses (COD, nitrogen species, phosphorus species and alkalinity) were conducted in accordance with the Standard Methods for the Examination of Water and Wastewater, APHA [23]. The pH of the reactors was determined using a DELTA 320 (Mettler Toledo, Columbus, OH, USA). Gas samples were taken from the reactors via a small sampling port located at the top of the reactor. The biogas production was measured with a wet gas meter (COMBIMASS @ GA-s hybrid premium, Binder Group, Germany, and the CH 4 production was determined using a gas chromatograph (Clarus 580 Arnel, PerkinElmer, Waltham, MA, USA) equipped with a thermal conductivity detector. Standard methods were used to measure TS, VS, TSS, and VSS [23]. Sludge samples were collected from the AnMBR and PSD and centrifuged at 10,000× g for 10 min at 4 • C. The remaining pellet of the sludge was transported on ice to the lab, where it was stored at −70 • C before DNA extraction.

DNA Extraction and 16S rRNA Gene Amplicon Analysis
We collected three samples from PSD and AnMBR, respectively. For each sample, three replicates of DNA extraction were conducted individually, and combined as a composition sample, before the PCR reaction. A Powersoil TM DNA Isolation Kit (Qiagen, Hilden, Germany) was used to extract microbial DNA from PSD and AnMBR according to the manufacturer's protocol. The yield and quality of DNA were analyzed electrophoretically on 1% agarose gel. Prokaryotic universal primers 515F and 806R [24], where the barcode was an 12-base sequence unique to each sample, were used in PCR to amplify the V4 region of the bacterial 16S ribosomal RNA gene. The PCRs were carried out in 50 µL volume with 10 ng of template DNA, 25 µL of 2 × ExTaq PCR Master Mix (Takara Biotechnology, Dalian, China) and 1 µL of each 10 mM primer. The PCR amplification was run under the following conditions: (1) pre-denature at 95 • C for 2 min; (2) 28 cycles of denaturing at 95 • C for 30 s, annealing at 55 • C for 30 s, and extension at 72 • C for 30 s; and (3) a final extension at 72 • C for 5 min. The PCR products were detected by 2% denaturing agarose gels, after which the 16S rRNA gene amplicons were extracted from the gels and purified with a QIAquick PCR purification kit (Qiagen, Düsseldorf, Germany) and quantified with Fluorometer (Qubit ® 3.0, ThermoFisher Scientific, Waltham, MA, USA). The purified 16S rRNA gene amplicons were paired-end sequenced on an Illumina Hiseq platform at Novogene Bioinformatics Technology, Beijing, China.
Analyses of the sequence of 16S rRNA genes were performed using QIIME2 software [25]. The following parameters were used for quality filtering: minimum/maximum length = 250/600; maximum 2 ambiguous bases; maximum 2 primer mismatches; average quality score > 25; and homopolymer < 6 nucleotides. OTUs with 97% sequence identity were clustered using the UCLUST software and the Greengenes, used as a reference 16S rRNA gene sequence database [26]. A representative sequence was chosen for each phylotype and aligned against the Greengenes using PyNAST [27]. The UCHIME algorithm using USEARCH (v7.0) was used to remove potentially chimeric sequences from the aligned representative sequences [28]. The taxonomic assignment of representative sequences from each OTUs was performed using the UCLUST consensus taxonomy assigner [29]. Finally, to estimate alpha diversity, a random sub-sampling method for each sequence library was used for microbial community diversity index calculations to control for the effects of library size. Alpha diversity indices were calculated for all samples with 1000 repetitions using a size of 16,126 sequences per sample. For principal coordinate analysis (PCoA), all samples were also subsampled to 16,126 sequences per sample to remove sample-size effects.

Metagenomic Sequencing, Assembly, and Annotation
The community genomic DNA (from triple DNA extractions from each replicate of the bioreactor) of the samples was sheared into 500-and 800-bp length fragments and multiplex-sequenced on the Illumina HiSeq 2500 platform (Illumina, Inc., San Diego, CA, USA), which generated 2 × 250-bp pair-end sequencing reads. Sequencing was carried out at Novogene Bioinformatics Technology (Beijing, China). The adapter sequences were removed from the generated reads, which were trimmed using Trimmomatic v0.30 with a quality cutoff of 30, a sliding window of 6 bp, and a minimum length cutoff of 45 bp [30]. The functional composition of the microbial communities was obtained from the results of the high-quality reads compared with the sequences in the Kyoto Encyclopedia of Genes and Genomes (KEGG, a database resource for understanding high-level functions and utilities of the biological system) [31] and Carbohydrate-Active EnZymes (CAZy, a database describes the families of structurally-related catalytic and carbohydrate-binding modules of enzymes) databases [32] using thresholds of e < e −10 and an aligned length > 150 nt. High-quality reads from the PSD and AnMBR were pooled and assembled into one metagenome using IDBA-UD 1.1.1 [33] with the following parameters: minimum k value of 60 and maximum k value of 120; increment of the k-mer of each iteration of 10 and minimum multiplicity for filtering k-mers when building the graph of 5; seed k-mer size for alignment of 5; and minimum length of contigs of 1000. The other parameters were set as default. The contigs and unmapped reads from the IDBA-UD assembly were re-assembled using the parameters as described. The contig coverage was determined by mapping the initial reads to contigs larger than 1 kbp with Bowtie2 software [34]. Subsequently, Samtools [35] was used to convert, sort, and merge the Sam files. Then, coverage (describes the number of reads that align to the contigs) was determined with the Genomecov software of the Bedtools package [36]. As low coverage and short contigs are known to be error-prone, contigs with a length < 1 kb and average coverage < 3 were discarded from the subsequent binning step.

Methanogenic Pathway of Predominant Species Determined by Genomic Binning
CONCOCT software was used for contig binning, which is based on sequence composition and coverage [37]. Subsequently, mapping of the high-quality reads was carried out using Bowtie2 [34] to determine the coverage of these contigs before conducting the CONCOCT step. CheckM [38] was used to assess the completeness and contamination of the recovered metagenome-assembled genomes (MAGs). MAGs with completeness > 90% and contamination < 10% were used for detailed downstream analyses. The abundance of each MAG of the samples was estimated by mapping the high-quality reads of individual datasets to the contigs from the MAGs.
Taxonomic analysis of the MAGs was performed by different methods, and the results were compared to obtain a consensus assignment. First, all the protein sequences translated in silico from these MAGs were input to PhyloPhlAn [39] to facilitate the phylogenetics of the MAGs, which were based on 400 universal proteins. Second, the ORFs predicted from each MAG were checked by sequence similarity to the NCBI-nr database using BlastX with a maximum allowed e-value of 1 × 10 −10 . The BlastX results were parsed by MEGAN [40] to access the species and genus phylum-level taxonomic annotation of each cluster of bins. The minimum bit score used for the analysis was 60, and a minimum support of 5% for each taxonomic category was used for the LCA algorithm. A taxon of a MAG was assigned when at least 85% of the identified ORFs resulted in a concordant taxonomy. We also constructed a phylogenetic tree of the recovered MAGs and their closest species' genome using the FastTree algorithm with default settings by PhyloPhlAn.
Genes were predicted across each MAG using FragGeneScan [41], and functional prediction of the identified genes was conducted by similarity searching against the KEGG and CAZy databases using an e-value of 10 −5 and an aligned length larger than 200 nt. The proteins of interest in carbohydrate and methanogenic metabolism were retrieved from the gene information in comparison with the KEGG database, and the abundance of methanogenesis-associated genes in each MAG was calculated using the mapping information of each gene.

Statistics
Statistical analyses were conducted using SPSS 17.0 Software (SPSS, Inc., Chicago, IL, USA). Statistical significance was estimated with two-way analysis of variance (ANOVA) followed by the LSD or Dunnett's T3 test to determine the significance of differences between groups. Linear discriminant analysis (LDA) of effect size (LEfSe) was applied to determine the most discriminant taxa among the PSD and AnMBR samples [42], with a value for the statistical test equal to 0.05 and a logarithmic LDA score threshold of 4.0.

Sequence Accession
The sequence data of the 16S rRNA genes and metagenome from this study are publicly available in NCBI BioProject PRJNA578561 (with accession no. SRR10566884-SRR10566893).

Operational Performance of the DWTS
The detailed performance of the DWTS is summarized in Table S1. The DWTS demonstrated good performance in terms of COD (52% on the average) and total N destruction (44%). The SO 4

2−
Appl. Sci. 2020, 10, 135 6 of 19 in the effluent of the system was significantly lower (2 mg/L) than that in the effluent (34.1 mg/L). A summary of the chemical characteristics of the PSD and AnMBR in the DWTS is provided in Table 1. There was a large difference in the organic loading rate between the AnMBR and PSD, in which 0.71 g of COD and 0.16 g of VSS per L wastewater were loaded into these two bioreactors every day. The DWTS functioned efficiently in regard to biogas production. Approximately 122 L per kg COD in the AnMBR and 374 L of biogas per kg volatile suspended solids (VSS) in the PSD were observed, and the average methane composition accounted for approximately 80%~85% and 60~65% of the biogas in the AnMBR and PSD, respectively. These results implied that a considerable amount of organic matter was converted into methane in both bioreactors. Similar pH, TSS, and VSS values were monitored in both reactors ( Table 1). The temperature was kept at approximately 20 • C and 35 • C for the AnMBR and PSD, respectively. In terms of methane and biogas recovery, the PSD offers potential advantages over the AnMBR. Table 1. Anaerobic membrane bioreactor (AnMBR) and sludge anaerobic reactor (PSD) system operating parameters.

Characteristic Microorganisms in the PSD and AnMBR
To identify the microorganisms present in the PSD and AnMBR, we conducted 16S rRNA gene amplicon sequencing on three samples from each bioreactor. After quality control, denoising and chimer removal, a total of 148,306 high-quality reads were generated from the 16S rRNA V4 sequences from 6 samples, ranging from 16,126 to 41,310 reads, with an average of 24,718 (Table S2). The sequence numbers were normalized to 16,126 reads in each sample to describe the alpha diversity and composition characteristics. All the sequences were clustered into OTUs at 97% sequence similarity. The number of OTUs per sample ranged between 386 and 608 (Table S2). Good's coverage ranged from 0.85595 to 0.93668 in all samples (Table S2), suggesting that major microbial diversity had been captured. The average Chao1 index values were 2217 ± 295 and 1459 ± 748 for the PSD and AnMBR samples, respectively. The PD indices were 57.4 ± 1.82 for the PSD samples and 49.5 ± 7.36 for the AnMBR samples. The Shannon and Simpson diversity indices were significantly higher in the AnMBR (Shannon: 7.012 ± 0.330; Simpson 0.978 ± 0.009) than in the PSD samples (Shannon: 5.973 ± 0.040; Simpson 0.942 ± 0.05) (p < 0.05) However, the number of observed OTUs, Chao1, ACE, and Faith PD index did not differ significantly between the PSD and AnMBR (Table S3).
The 16S rRNA gene sequences were classified at both the phylum and family levels. At the phylum level, Proteobacteria (27.11 ± 0.46%), Euryarchaeota (13.72 ± 0.41%), Firmicutes (11.44 ± 0.92%), Synergistetes (11.22 ± 1.02%), and Bacteroidetes (10.89 ± 2.51%) were the predominant phyla in the PSD samples, accounting for 74.38% of the total microbiota. However, a large proportion of Proteobacteria (75.53 ± 0.46%), Firmicutes (6.37 ± 0.92%), and Bacteroidetes (5.76 ± 2.51%) were observed in the AnMBR samples ( Figure 2). Samples of the PSD and AnMBR also harbored different family-level profiles, with Methanobacteriaceae (12.99 ± 0.35% in the PSD and 2.08 ± 0.35% in the AnMBR) and Syntrophaceae (4.24 ± 0.10% in the PSD and 3.61 ± 0.01% in the AnMBR) representing the most abundant genera in the PSD and AnMBR samples, respectively (Table S4). At the phylum level, the relative abundances of most predominant phyla were significantly different between the two types of samples (all p < 0.05) ( Figure 2). Moreover, Methanobacteriaceae (12.99 ± 0.35% in the PSD and 2.08 ± 0.35% in the AnMBR) and Syntrophaceae (4.24 ± 0.10% in the PSD and 3.61 ± 0.01% in the AnMBR) were identified as the families that most significantly differed in abundance between the PSD and AnMBR samples (all p < 0.05) (Table S4). This result suggested a distinct difference in the microbial community structure between the PSD and AnMBR samples, which was also supported by the weighted PCoA results. The weighted PCoA showed that the PSD samples grouped together in one cluster, while the AnMBR samples grouped in another cluster ( Figure S1). Methanobacteriaceae and Syntrophaceae were identified as the genera that significantly differed in abundance between the PSD and AnMBR samples. Methanobacteriaceae is recognized as a family comprising hydrogenotrophic methanogenic archaeon in anaerobic digestion and is prominent in methanogenic communities [43]. The hydrogenotrophic Methanobacteriales have been shown to correlate with biogas production in 29 full-scale digester studies, confirming their role in the maintenance of digester function [44]. Methanobacteria use the methanogenesis pathway, in which hydrogen reacts with carbon dioxide to produce methane [45]. The maintenance of the temperature (PSD at 37 • C and AnMBR at 20 • C) should be the primary reason for the significantly different abundance of Methanobacteriaceae in these two bioreactors and should also explain why biogas production was higher in the PSD. Compared with the PSD, the relative abundance of Syntrophaceae was significantly higher in the AnMBR. Bacteria affiliated with Syntrophaceae play important roles in acetate oxidation during anaerobic methane production and are seen as stabilizers that maintain AD stability, especially when the system undergoes environmental fluctuations [46]. Furthermore, the significance of Syntrophaceae in the methanogenic degradation of waste has been acknowledged [47,48]. The remarkably different abundances of Syntrophaceae and Methanobacteriaceae may reflect the functional tendencies in the AnMBR and PSD, where acetate oxidation and methane generation are the dominant microbial metabolisms, respectively, in these two types of bioreactors during wastewater treatment. The relative abundances of Desulfomonile, Smithella, and Syntrophus (all within the family Syntrophaceae) were significantly different between samples of PSD and AnMBR (p < 0.01). However, the above three genera were less abundant (<1%) in these two kinds of bioreactors. In the family Methanobacteraceae, Methanobacterium was the only identified genus with an averagely relative abundance larger than 1% in PSD and AnMBR, and also, the identified genus exhibited significantly different abundance in these two kinds of bioreactors.
relative abundances of Desulfomonile, Smithella, and Syntrophus (all within the family Syntrophaceae) were significantly different between samples of PSD and AnMBR (p < 0.01). However, the above three genera were less abundant (<1%) in these two kinds of bioreactors. In the family Methanobacteraceae, Methanobacterium was the only identified genus with an averagely relative abundance larger than 1% in PSD and AnMBR, and also, the identified genus exhibited significantly different abundance in these two kinds of bioreactors. Taxonomic composition of the microbial communities showing the differences between the PSD and AnMBR samples. The data reveal the phylum-and genus-level classification of the 16S rRNA gene sequences. 16S rRNA gene sequences longer than 250 bp were classified using the RDP naïve Bayesian rRNA Classifier at an 80% confidence threshold. Only the relative abundances of the identified phyla higher than 1% are listed. The average ± SD values of the samples in each group are expressed in each column. The significance of the cutoff for the corrected p-value using the Benjamini-Hochberg FDR procedure was set at 0.05. The significant differences within the sampling group are labeled **, where p < 0.05.

Variation of Gene Functional Profiles of the PSD and AnMBR
Based on the above results regarding the geochemical properties and microbial communities of the PSD and AnMBR, we hypothesize that the microbial interactions and metabolic pathways of the dominant species are very different in the two bioreactors. In addition, most of the 16S rRNA gene sequences could not be classified at the family level, indicating most of the species in the PSD and AnMBR are less well-understood. Thus, comparative metagenomics was carried out to explore the divergence of the predicted microbial metabolism during methane generation in these two bioreactors. Additionally, MAGs were reconstituted to reveal the interactions of the syntrophic consortium during methane generation. Large data sets of metagenomic DNA were produced for the PSD and AnMBR samples. We generated a total of 35.9 million 2 × 250 bp paired-end reads, with 10.34 G and 7.27 G nucleotide bases retained after quality filtering for the PSD and AnMBR samples, respectively.
The metabolic potential of the microbial communities of the PSD and AnMBR was predicted by comparison to the KEGG Ortholog (KO) database. In the PSD and AnMBR metagenomes, genes encoding metabolic functions were dominant, representing 20.22% and 22.79% within each sample, respectively, followed by genetic information processing (4.56% in the PSD and 5.07% in the AnMBR) and environmental information processing (3.30% in the PSD and 4.25% in the AnMBR) ( Figure S2). In the category of metabolism, the most abundant metabolic type was carbohydrate metabolism (4.28% in the PSD and 4.82% in the AnMBR), followed by amino acid metabolism (4.19% in the PSD and 4.73% in the AnMBR) as well as energy metabolism (3.15% in the PSD and 3.35% in the AnMBR) ( Figure S2). A high proportion of genes related to these functions has been previously reported in metagenomes [49][50][51], metatranscriptomes [52] and metaproteomes [53] from AD systems, indicating these metabolic activities are linked to the conversion of excess carbohydrates into methane during anaerobic digestion. When classified as KEGG level 3, the abundance of genes in methane metabolism was one of the categories with the highest abundance in the PSD and AnMBR, consisting of 1.05% Appl. Sci. 2020, 10, 135 9 of 19 and 0.88% of the total sequences (Table S5), respectively, suggesting methanogenesis was the major microbial metabolism in these two bioreactors.
To better understand the difference in carbohydrate decomposition in the PSD and AnMBR, we screened our metagenome for reads annotated as CAZyme. The PSD metagenome contained 864,945 reads (consisted of 4.17% of the total reads) associated with CAZyme that were unequally distributed among glycosyl transferases (GTs), glycoside hydrolases (GHs), carbohydrate binding modules (CBMs), carbohydrate esterases (CEs), auxiliary activities (AAs), and polysaccharide lyases (PLs) (Table S6), and their read abundances were 41.3%, 37.2%, 12.8%, 5.3%, 2.6%, and 0.8%, respectively ( Figure S3). In the AnMBR metagenome, the read abundances of GTs, GHs, CBMs, CEs, AAs, and PLs were 43.4%, 37.8%, 12.0%, 4.4%, 2.6%, and 0.8%, respectively ( Figure S3). A higher gene abundance of CAZyme suggests a greater role of carbohydrate metabolism in the anaerobic digestion of the PSD and AnMBR. GTs and GHs had higher percentages in terms of total reads, indicating glycosyl transferases and glycoside hydrolases play key roles in the cleavage of polymeric substrates in DWTS. The PSD metagenomes recorded a diverse profile of carbohydrate pathways for a total of 252 GHs, 85 GTs, 67 PLs, 16 CEs, 64 CBMs, and 20 AAs families, comparable to 254, 88, 66, 16, 69, and 21 in the AnMBR metagenome, respectively (Table S6). Approximately 357,389 and 321,643 reads were assigned to different GT and GH families in the PSD compared with 304,909 GTs and 258,444 GHs in the AnMBR metagenome (Table S6), implying their common abilities to form and hydrolyze glycosidic bonds when breaking down polysaccharides in the DWTS [54].

Phylogenetics of MAGs and Their Abundance in the PSD and AnMBR
A large fraction of the reads were assembled in contigs ≥1000 bp (114,558 contigs including 396,835,232 bp), yielding an N50 (defined as the minimum contig length needed to cover 50% of the total lengths of cotnigs) of 4676 with a max contig length of 323,550 and a mean size of contigs of 2139 bp. To further examine the metabolic properties of the dominant species in the PSD and AnMBR, metagenomic binning based on composition and differential coverage of contigs was performed to reconstruct the MAGs of the species that dominated these processes. Thereby, 21 MAGs (16 and 7 affiliated with Bacteria and Archaea, respectively) were observed and assigned to 15 categories, with an average completeness of 93.0% and contamination of 2.9% ( Table 2). The genome sizes of the MAGs ranged from 1,551,831 bp to 5,015,715 bp, with an average of 2,650,305 bp, and the number of contigs in each MAG ranged from 142 to 629, with an average N50 size of 12,914 and mean contig length of 8695 bp (Table 2). A detailed overview of the taxonomic classification of the MAGs is provided in Figure 3 and Table 2, in which a taxonomic classification with high confidence was possible to the species level for all the MAGs. In summary, 7 MAGs were attributed to the phylum Euryarchaeota,  (Figure 3). The relationships between the MAGs and their closest species with complete genome sequences are shown in Figure 3. Most of the recruited MAGs belonged to Methanogenic archaea. MAG13, MAG14, MAG16, and MAG19 were most familiar with the genome of Methanobacterium sp. SMA-27, for which no information regarding its origin is publicly available. The sequences of MAG215 had high similarity with Methanomassiliicoccales archaeon RumEn M2, which has been reported to be a trimethylamine-utilizing archaea [55]. The closest genome affiliated with MAG11 was Methanospirillum hungatei JF-1, whereas the genome of MAG220 was affiliated with Methanosaeta harundinacea 6Ac. During the degradation of complex organic matter to carbon dioxide and methane in anaerobic environments, M. hungatei uses hydrogen or formate to produce methane, which are produced by other syntrophic microbes [56]. M. harundinacea is an acetate-scavenging methanogen isolated from a upflow anaerobic sludge blanket reactor [57]. A Bacteroidetes endosymbiont of Geopemphigus sp. exhibited the highest relationship with MAG158, MAG217, MAG228, MAG42, and MAG83. Prolixibacteraceae bacterium XSD2, Fusobacteriaceae bacterium, Nitrospirae bacterium HCH-1, and Desulfovibrio cuneatus DSM 11,391 were most similar with the genomes of MAG59, MAG142, MAG161, and MAG68, respectively. D. cuneatus is a psychrotolerant sulfate-reducing bacterium isolated from an oxic freshwater sediment [58]. In addition, MAG55, MAG100, MAG74, MAG200, and MAG159 showed the highest similarity with Sulfurovum sp. NBC37-1, Succinatimonas hippei YIT 12066, Aminiphilus circumscriptus DSM 16,581, Synergistaceae bacterium, and Succinispira mobilis DSM 6222, respectively. Sulfurovum sp. NBC37-1 is a chemolithoautotroph deep-sea Epsilonproteobacteria that obtains energy through mesophilic hydrogen-and sulfur-oxidation [59]. S. hippei YIT 12066T was isolated from human feces, and the growth of S. hippei YIT 12,066 T depends on CO 2 or bicarbonate [60]. A. circumscriptus is an anaerobic amino-acid-degrading bacterium from an upflow anaerobic sludge reactor [61], and S. mobilis was proven to be a succinate-decarboxylating anaerobic bacterium in a mixed culture growing with glycolate [62]. To estimate the relative abundance of the MAGs represented in the PSD and AnMBR, all the quality-filtered reads from each sample were aligned to the contigs from each MAG with Bowtie2 software. As shown in Figure 3, the relative abundance of the MAGs exhibited a large divergence between the PSD and AnMBR. For example, the highest coverage was observed in MAG16 (affiliated with Methanobacterium sp., with 5.27% reads mapped on the corresponding contigs in PSD), followed by MAG200 (affiliated with Synergistaceae bacterium, 2.79%), MAG13 (affiliated with Methanobacterium sp.), MAG74 (affiliated with A. circumscriptus) and MAG220 (affiliated with M. concilii), which were the other predominant species in the PSD, with abundances of 1.42%, 1.18%, and 0.96%, respectively. In the AnMBR metagenome, MAG83 (classified as Bacteroidetes bacterium), MAG200 (Synergistaceae bacterium), MAG217 (Bacteroidetes bacterium), MAG42 (Bacteroidales bacterium), and MAG59 (Prolixibacteraceae bacterium) were the most abundant species, with relative abundances of 4.74%, 2.93%, 2.85%, 1.59%, and 1.36%, respectively ( Figure 3). Methanobacteriaceae is capable of directly converting gaseous hydrogen and carbon dioxide into methane [63], whereas microorganisms belonging to Synergistaceae family (Aminobacterium thunnarium, Aminobacterium colombiense, and Aminivibrio pyruvatiphilus) are certified as syntrophic acetate oxidizers [64,65]. This result may also derive from the functional difference between the PSD and AnMBR in the DWTS, in which wastewater was solid-liquid separated, and organic matter was degraded in the AnMBR, followed by methanogenesis in the sludge of PSD bioreactors.

Microbial Interaction and Methane-Producing Pathways of the Dominant Species in the PSD and AnMBR
To better understand the predicted function carried out by each MAG of the draft genome, genes were annotated through BlastX analyses against the KEGG and CAZy databases. The number of predicted ORFs in MAGs ranged from 1658 to 4635 with an average of 2648 (Table 2). On average, 50.8% (ranged from 36.4% to 63.7%) and 2.44% (ranged from 0.74% to 8.66%) of the ORFs were annotated in the KEGG and CAZy databases (Table S7). Genes associated with carbohydrate-active and methanogenic metabolism were recruited from each MAG, benefiting the prediction of the functions carried out by the predominant species during methane generation. The abundance of CAZyme-encoding genes varied across the MAGs. MAG59 (classified as Prolixibacteraceae bacterium), MAG158 (Bacteroidetes bacterium) and MAG127 (Bacteroidetes bacterium) had the most To estimate the relative abundance of the MAGs represented in the PSD and AnMBR, all the quality-filtered reads from each sample were aligned to the contigs from each MAG with Bowtie2 software. As shown in Figure 3, the relative abundance of the MAGs exhibited a large divergence between the PSD and AnMBR. For example, the highest coverage was observed in MAG16 (affiliated with Methanobacterium sp., with 5.27% reads mapped on the corresponding contigs in PSD), followed by MAG200 (affiliated with Synergistaceae bacterium, 2.79%), MAG13 (affiliated with Methanobacterium sp.), MAG74 (affiliated with A. circumscriptus) and MAG220 (affiliated with M. concilii), which were the other predominant species in the PSD, with abundances of 1.42%, 1.18%, and 0.96%, respectively. In the AnMBR metagenome, MAG83 (classified as Bacteroidetes bacterium), MAG200 (Synergistaceae bacterium), MAG217 (Bacteroidetes bacterium), MAG42 (Bacteroidales bacterium), and MAG59 (Prolixibacteraceae bacterium) were the most abundant species, with relative abundances of 4.74%, 2.93%, 2.85%, 1.59%, and 1.36%, respectively ( Figure 3). Methanobacteriaceae is capable of directly converting gaseous hydrogen and carbon dioxide into methane [63], whereas microorganisms belonging to Synergistaceae family (Aminobacterium thunnarium, Aminobacterium colombiense, and Aminivibrio pyruvatiphilus) are certified as syntrophic acetate oxidizers [64,65]. This result may also derive from the functional difference between the PSD and AnMBR in the DWTS, in which wastewater was solid-liquid separated, and organic matter was degraded in the AnMBR, followed by methanogenesis in the sludge of PSD bioreactors.

Microbial Interaction and Methane-Producing Pathways of the Dominant Species in the PSD and AnMBR
To better understand the predicted function carried out by each MAG of the draft genome, genes were annotated through BlastX analyses against the KEGG and CAZy databases. The number of predicted ORFs in MAGs ranged from 1658 to 4635 with an average of 2648 (Table 2). On average, 50.8% (ranged from 36.4% to 63.7%) and 2.44% (ranged from 0.74% to 8.66%) of the ORFs were annotated in the KEGG and CAZy databases (Table S7). Genes associated with carbohydrate-active and methanogenic metabolism were recruited from each MAG, benefiting the prediction of the functions carried out by the predominant species during methane generation. The abundance of CAZyme-encoding genes varied across the MAGs. MAG59 (classified as Prolixibacteraceae bacterium), MAG158 (Bacteroidetes bacterium) and MAG127 (Bacteroidetes bacterium) had the most abundant CAZy genes, with total numbers of 269, 232, and 103, respectively. In carbohydrate-active modules, MAG59 had the highest abundance of genes encoding GHs, GTs, and PLs, with gene numbers of 191, 46, and 10 ( Figure 4), respectively, providing genetic evidence for the high cellulolytic and transglycosylatic abilities of Prolixibacteraceae bacterium. The most abundant genes of CEs and CBMs were identified in MAG158. The above results suggest that polysaccharides within the DWTS habitat are hydrolyzed by the predominant Bacteroidetes bacterium and Prolixibacteraceae bacterium.
Appl. Sci. 2020, 10, x FOR PEER REVIEW abundant CAZy genes, with total numbers of 269, 232, and 103, respectively. In carbohydrate-active modules, MAG59 had the highest abundance of genes encoding GHs, GTs, and PLs, with gene numbers of 191, 46, and 10 ( Figure 4), respectively, providing genetic evidence for the high cellulolytic and transglycosylatic abilities of Prolixibacteraceae bacterium. The most abundant genes of CEs and CBMs were identified in MAG158. The above results suggest that polysaccharides within the DWTS habitat are hydrolyzed by the predominant Bacteroidetes bacterium and Prolixibacteraceae bacterium.  Figure S5). For hydrogenotrophic methanogenesis, carbon dioxide is successively reduced to methane through a series of intermediates with a methyl group. The methyl group is then transferred to Coenzyme M, forming methyl-CoM. Thus, methyl-CoM is reduced to methane through methyl coenzyme M reductase (Mcr) at the final step. Genes encoding enzymes that participate in hydrogenotrophic methanogenesis had the highest abundance in the archaeal MAGs, followed by methylotrophic and acetotrophic methanogenesis ( Figure 5, Figures S4 and S5). Generally, the genes involved in the hydrogenotrophic methanogenesis pathway had higher abundance in archaea MAGs ( Figure 5), especially for formylmethanofuran dehydrogenase (K00200-K00203) and tetrahydromethanopterin S-methyltransferase (K00577-K00584). The genes of MAG14 and MAG16 (both classified as Methanobacterium sp.) included almost all the enzymes in the completely hydrogenotrophic pathway, indicating Methanobacterium-associated archaea tended to  Figure S5). For hydrogenotrophic methanogenesis, carbon dioxide is successively reduced to methane through a series of intermediates with a methyl group. The methyl group is then transferred to Coenzyme M, forming methyl-CoM. Thus, methyl-CoM is reduced to methane through methyl coenzyme M reductase (Mcr) at the final step. Genes encoding enzymes that participate in hydrogenotrophic methanogenesis had the highest abundance in the archaeal MAGs, followed by methylotrophic and acetotrophic methanogenesis ( Figure 5, Figures S4 and S5). Generally, the genes involved in the hydrogenotrophic methanogenesis pathway had higher abundance in archaea MAGs ( Figure 5), especially for formylmethanofuran dehydrogenase (K00200-K00203) and tetrahydromethanopterin S-methyltransferase (K00577-K00584). The genes of MAG14 and MAG16 (both classified as Methanobacterium sp.) included almost all the enzymes in the completely hydrogenotrophic pathway, indicating Methanobacterium-associated archaea tended to convert carbon dioxide to methane using hydrogen, which has been proven by genomic and physiological data [68,69]. It was interesting that the genes for complete typical hydrogenotrophic pathway were not detected in MAG215 (classified as Methanomassiliicoccales archaeon), while MAG215 harbored a high abundance of genes involved in methylotrophic methanogenesis; both of these results suggest that Methanomassiliicoccales archaea produced methane from methanol, in agreement with genomic insights into the methylotrophic lifestyle of Methanomassiliicoccales [70]. The acetotrophic pathway is known to be the major pathway for more than 70% of methane production in most engineering anaerobic digestion processes [71]. In the acetotrophic methanogenesis pathway ( Figure S4), genes encoding acetyl-CoA synthetase (K01895) and acetyl-CoA decarbonylase/synthase (K00193, K00194, K00197) are present in all archaeal MAGs and are involved in acetate oxidation [72]. However, genes encoding phosphate acetyltransferase (K00625 and K13788) and acetate kinase (K00925), which transform CoA into acetyl-CoA and acetate into acetyl phosphate, respectively, were absent. Furthermore, no archaeal MAGs harbored all the genes encoding enzymes in the complete pathway of acetotrophic methanogenesis, suggesting the predominant pathway for methane generation was not through acetate conversion in the PSD and AnMBR bioreactors. In the methylotrophic pathway ( Figure S5), almost all the archaeal MAGs (except MAG11, classified as Methanomicrobiales archaeon) had genes encoding the three subunits of methyl-coenzyme M reductase, which is the key enzyme of biological methane formation [73]. However, all the subunits of heterodisulfide reductase have not been identified in all of these archaea MAGs, posing a question regarding the ability of the affiliated archaea to convert methanol to methane. However, it should be noted that the methanogenesis pathway was predicted based on the abundance of genes in MAGs, rather than metatranscriptomics or metaproteomics, which are required to further explore the active functions involved in the methanogenesis pathway in future studies.
MAG215 harbored a high abundance of genes involved in methylotrophic methanogenesis; both of these results suggest that Methanomassiliicoccales archaea produced methane from methanol, in agreement with genomic insights into the methylotrophic lifestyle of Methanomassiliicoccales [70]. The acetotrophic pathway is known to be the major pathway for more than 70% of methane production in most engineering anaerobic digestion processes [71]. In the acetotrophic methanogenesis pathway ( Figure S4), genes encoding acetyl-CoA synthetase (K01895) and acetyl-CoA decarbonylase/synthase (K00193, K00194, K00197) are present in all archaeal MAGs and are involved in acetate oxidation [72]. However, genes encoding phosphate acetyltransferase (K00625 and K13788) and acetate kinase (K00925), which transform CoA into acetyl-CoA and acetate into acetyl phosphate, respectively, were absent. Furthermore, no archaeal MAGs harbored all the genes encoding enzymes in the complete pathway of acetotrophic methanogenesis, suggesting the predominant pathway for methane generation was not through acetate conversion in the PSD and AnMBR bioreactors. In the methylotrophic pathway ( Figure S5), almost all the archaeal MAGs (except MAG11, classified as Methanomicrobiales archaeon) had genes encoding the three subunits of methyl-coenzyme M reductase, which is the key enzyme of biological methane formation [73]. However, all the subunits of heterodisulfide reductase have not been identified in all of these archaea MAGs, posing a question regarding the ability of the affiliated archaea to convert methanol to methane. However, it should be noted that the methanogenesis pathway was predicted based on the abundance of genes in MAGs, rather than metatranscriptomics or metaproteomics, which are required to further explore the active functions involved in the methanogenesis pathway in future studies. Various microorganisms are involved in and cooperate with each other to achieve methane formation in anaerobic digestion. According to the genome binning presented in this work, members of Bacteroidetes and Prolixibacteraceae were predicted to hydrolyze polysaccharide to acetate. Subsequently, Synergistaceae may oxidize acetate to hydrogen and carbon dioxide. Finally, Methanobacterium sp. was possible to reduced carbon dioxide to methane with hydrogen as the primary electron donor. Furthermore, sulfate was likely to be removed through desulfurization with the help of Desulfovibrio desulfuricans ( Figure 6).
Appl. Sci. 2020, 10, x FOR PEER REVIEW Various microorganisms are involved in and cooperate with each other to achieve methane formation in anaerobic digestion. According to the genome binning presented in this work, members of Bacteroidetes and Prolixibacteraceae were predicted to hydrolyze polysaccharide to acetate. Subsequently, Synergistaceae may oxidize acetate to hydrogen and carbon dioxide. Finally, Methanobacterium sp. was possible to reduced carbon dioxide to methane with hydrogen as the primary electron donor. Furthermore, sulfate was likely to be removed through desulfurization with the help of Desulfovibrio desulfuricans ( Figure 6).

Conclusions
The in-depth 16S rRNA gene amplicon and metagenome information obtained in this study contribute to the understanding of the taxonomies and pathways involved in methanogenesis in the two types of bioreactors in the DWTS. The results obtained in this study showed that microbial composition was significantly different between the PSD and AnMBR. High-quality MAGs of predominated species were recruited, and the results of an in-depth dissection of the genes in the MAGs suggested a strong interaction of bacteria and archaea during methanogenesis that mainly occurs through the interspecific exchange of substances. The valuable insights obtained in this study are expected to provide guidance for the currently planned start-up of decentralized wastewater treatment systems that contain a primary sludge digestor and membrane bioreactor.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1. Table S1, physicochemical characteristics of the influent and effluent water; Table S2, sequence numbers and alpha diversity indices of 16S rRNA genes in the PSD and AnMBR samples; Table S3, differences of the alpha diversity indices of the 16S rRNA genes between the PSD and AnMBR samples; Table S4, taxonomic profiling of the 16S rRNA genes at the family level in the PSD and AnMBR; Table S5, KEGG categories of the PSD and AnMBR metagenomes; Table S6, number of reads annotated as genes affiliated with CAZyme in the PSD and AnMBR metagenomes; Table S7, rate of ORFs in MAGs annotated by the KEGG and CAZy databases; Figure  S1, weighted PCoA clustering of the microbial communities in the PSD (blue triangles) and AnMBR (red triangles) samples based on the UniFrac distance; Figure S2, relative abundances of the KEGG categories of functional genes in the PSD and AnMBR metagenomes; Figure S3, read abundance of CAZyme in the PSD and AnMBR metagenomes; Figure S4, numbers of genes involved in the relevant acetotrophic methanogenesis pathway in archaeal MAGs; Figure S5, numbers of genes involved in the relevant methylotrophic methanogenesis pathway in archaeal MAGs; Figure S6, the relative abundances of Syntrophaceae and Methanobacteriaceae in PSD and AnMBR.

Conclusions
The in-depth 16S rRNA gene amplicon and metagenome information obtained in this study contribute to the understanding of the taxonomies and pathways involved in methanogenesis in the two types of bioreactors in the DWTS. The results obtained in this study showed that microbial composition was significantly different between the PSD and AnMBR. High-quality MAGs of predominated species were recruited, and the results of an in-depth dissection of the genes in the MAGs suggested a strong interaction of bacteria and archaea during methanogenesis that mainly occurs through the interspecific exchange of substances. The valuable insights obtained in this study are expected to provide guidance for the currently planned start-up of decentralized wastewater treatment systems that contain a primary sludge digestor and membrane bioreactor.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/10/1/135/s1. Table S1, physicochemical characteristics of the influent and effluent water; Table S2, sequence numbers and alpha diversity indices of 16S rRNA genes in the PSD and AnMBR samples; Table S3, differences of the alpha diversity indices of the 16S rRNA genes between the PSD and AnMBR samples; Table S4, taxonomic profiling of the 16S rRNA genes at the family level in the PSD and AnMBR; Table S5, KEGG categories of the PSD and AnMBR metagenomes; Table S6, number of reads annotated as genes affiliated with CAZyme in the PSD and AnMBR metagenomes; Table S7, rate of ORFs in MAGs annotated by the KEGG and CAZy databases; Figure S1, weighted PCoA clustering of the microbial communities in the PSD (blue triangles) and AnMBR (red triangles) samples based on the UniFrac distance; Figure S2, relative abundances of the KEGG categories of functional genes in the PSD and AnMBR metagenomes; Figure S3, read abundance of CAZyme in the PSD and AnMBR metagenomes; Figure S4, numbers of genes involved in the relevant acetotrophic methanogenesis pathway in archaeal MAGs; Figure S5, numbers of genes involved in the relevant methylotrophic methanogenesis pathway in archaeal MAGs; Figure S6, the relative abundances of Syntrophaceae and Methanobacteriaceae in PSD and AnMBR.