Biotic and abiotic factors acting on community assembly in parallel anaerobic digestion systems from a brewery wastewater treatment plant

ABSTRACT

Anaerobic digestion is a complex microbial process that mediates the transformation of organic waste into biogas. The performance and stability of anaerobic digesters relies on the structure and function of the microbial community. In this study, we asked whether the deterministic effect of wastewater composition outweighs the effect of reactor configuration on the structure and dynamics of anaerobic digester archaeal and bacterial communities. Biotic and abiotic factors acting on microbial community assembly in two parallel anaerobic digestion systems, an upflow anaerobic sludge blanket digestor (UASB) and a closed digester tank with a solid recycling system (CDSR), from a brewery WWTP were analysed utilizing 16S rDNA and mcrA amplicon sequencing and genome-centric metagenomics. This study confirmed the deterministic effect of the wastewater composition on bacterial community structure, while the archaeal community composition resulted better explained by organic loading rate (ORL) and volatile free acids (VFA). According to the functions assigned to the differentially abundant metagenome-assembled genomes (MAGs) between reactors, CDSR was enriched in genes related to methanol and methylamines methanogenesis, protein degradation, and sulphate and alcohol utilization. Conversely, the UASB reactor was enriched in genes associated with carbohydrate and lipid degradation, as well as amino acid, fatty acid, and propionate fermentation. By comparing interactions derived from the co-occurrence network with predicted metabolic interactions of the prokaryotic communities in both anaerobic digesters, we conclude that the overall community structure is mainly determined by habitat filtering.

GRAPHICAL ABSTRACT

KEYWORDS:

• Community assembly
• microbial interactions
• anaerobic digestion
• metabolic modelling
• metagenome

Introduction

Anaerobic digestion (AD) is the method of choice for the treatment of high-strength wastewater. The basis of this technology involves organic matter degradation achieved by the concerted action of several metabolically interconnected microbial guilds. Bacteria carry out the initial stages of hydrolysis, acidogenesis (fermentation), and acetogenesis, whereas the final step, methanogenesis, is performed by methanogenic archaea [1]. It has been proposed that the performance and stability of anaerobic digesters are highly dependent on the microbial community structure. Information regarding the microbial community composition in AD has been increasing in recent years. However, there is still a lack of understanding respecting how to harness environmental variables to effectively manage these microbial communities with the aim of improving resilience, stability, and recovery after process alterations. For instance, performance differences in AD have been related to variations in microbial community composition [2]. In addition, efficiency has been associated with diversity [3], whereas functional stability has been related to richness and evenness [1,3,4]. Thus, microbial communities with higher functional diversity and higher evenness are presumed to exhibit greater resistance to external disturbances, such as changes in wastewater composition or operational parameters.

Several studies of wastewater treatment plants (WWTP) have highlighted the importance of deterministic factors in community assembly [5–7]. Specifically, wastewater composition has been pointed out as a key factor determining bacterial composition in WWTPs [8,9]. Ibarbalz et al. have demonstrated that despite differences in plant configuration and operation, WWTPs treating similar wastewater in different parts of the world, exhibit similar bacterial composition and genetic diversity [8]. Anaerobic microbiomes also appeared to be defined by operational conditions [10,11], including substrate characteristics [9,12]⁠. Similarly, reactor configuration has been proposed as a determinant of archaeal composition in AD systems [13]. Moreover, other deterministic factors, such as microbial interactions, have been proposed as major drivers of microbial community composition [14]. For example, it has been shown that increased competition among bacteria could explain co-occurrence patterns [15], whereas metabolic dependencies were identified as major drivers of species co-occurrence in simulated high-order communities [16].

In this study, we asked whether the deterministic effect of wastewater composition outweighs the effect of reactor configuration on the structure and dynamics of anaerobic digester archaeal and bacterial communities. To achieve this objective, we assessed the microbial community structure and dynamics of two parallel anaerobic reactors differing in configuration, treating wastewater from a full-scale industrial WWTP. The two reactors consisted of an Upflow Anaerobic Sludge Blanket (UASB) and a closed digester tank with a solid recycle system (CDSR) that received the same wastewater from a brewery during a two year period. We analysed the composition and dynamics of the prokaryotic communities present in both anaerobic digestion systems using amplicon sequencing targeting 16S rDNA for bacteria and mcrA gene, encoding the alpha subunit of the methyl-coenzyme M reductase, for methanogenic archaea. Firstly, we examined the relationship between microbial community structure and abiotic operational parameters. Secondly, we computed co-occurrence networks from taxa abundance data to infer microbial interactions participating as the most important biotic forces in both systems. Lastly, we employed metagenomic-based metabolic modelling to compare the effect of microbial interactions and habitat filtering on the assembly of their respective communities.

Our hypothesis was that despite the differences in reactor configuration, both reactors would share a similar bacterial and archaeal composition. Furthermore, we anticipated that there would be a relationship between operational conditions, microbial structure, stability, and efficiency of the wastewater treatment system.

Materials and methods

Sampling and data collection

Samples were collected from a brewery wastewater treatment plant located in the province of Buenos Aires, Argentina. During treatment, the raw effluent is first received by two equalizers for homogenization, before being split into the two anaerobic reactors, where most organic matter is degraded. The two anaerobic digesters have different scales and configurations. They operate in parallel during the peak production season, while one of them receives little or no wastewater during the winter period when beer production is low.

The largest reactor is a closed digester tank with a solid recycling system (CDSR), with a total volume of 1280 m³, that operates continuously year-round. It is compartmentalized into two sectors: the contact tank, where the wastewater comes into contact with the microorganisms, and a settling tank, where the solids settle and the sludge is returned to the contact tank in the downstream zone.

The other reactor, with a volume of 400 m³, is of the UASB (Upflow Anaerobic Sludge Blanket) type. In this type of configuration, wastewater enters the tank from the bottom and follows an upward flow thanks to the pressure exerted by the inlet flow. In the upper region of the tank, there are gas baffles and a solid–liquid separator so that the three phases (liquid effluent, solid sludge blanket and biogas) are separated.

Samples of anaerobic sludge were taken from both CDSR and UASB reactors on the dates indicated in Supplemental Figure 1, resulting in a total of 6 sampling instances over two years. In certain instances, more than one sample was taken per reactor to evaluate the variability within the reactors. Samples were transferred at room temperature within 2 h to the laboratory, where they were divided into 2-ml aliquots and stored at −20°C.

Operating parameter data were obtained from daily measurements at the WWTP conducted by plant personnel using standard methods. The analysed variables were pH, temperature, acidity, alkalinity, chemical oxygen demand (COD). Reactor performance efficiency was measured as the percent of COD removal. Volatile fatty acids (VFA) concentration was estimated from intermediate alkalinity [17,18]. For multivariate analysis, operational data were averaged from the values of the week before each sampling date.

DNA extraction and sequencing

Metagenomic DNA extraction was performed starting from a 500 μl aliquot of the homogenized sample using the FastDNA Spin Kit for Soil (MP Biomedicals, USA), according to the manufacturer’s instructions. V3-V4 16S rDNA amplicon sequencing was performed on an Illumina MiSeq platform at Macrogen Inc. with primers b341F (CCTACGGGAGGCAGCAG) and Bakt805R (GACTACHVGGGTATCTAATCC). For mcrA sequencing, samples were amplified in-house with primers Mlf (GGTGGTGTMGGATTCACACARTAYGCWACAGC) and Mlr (TTCATTGCRTAGTTWGGRTAGTT) [19], and sequenced at Macrogen Inc. with Illumina MiSeq tecnology. For shotgun sequencing, NextEra libraries of 937 bp on average were run as PE 2 × 250 pb in the HiSeq platform at INDEAR (Rosario, Argentina).

Raw reads from 16S rRNA and mcrA genes obtained in this study, and unassembled reads from the shotgun sequencing have been deposited in the Sequence Read Archive from the National Center for Biotechnology Information (NCBI; BioProject ID: PRJNA947194; SRR23930591/622 and SRR23949945/8 respectively).

Amplicon sequence processing

For amplicon sequencing, primers, poor quality reads, and short reads were trimmed with Trimmomatic [20]. Sequences were processed with the USEARCH v9.2.64 pipeline [21]. Chimeric sequences and those with expected error ≥1 were discarded. The remaining reads were mapped to OTUs with a minimum of 97% identity for 16SrDNA reads or 84% for mcrA reads [22] and the OTU tables were constructed. Data were rarefied to 14,250 and 27,224 sequences by sample for 16S rDNA and mcrA, respectively. Bacterial OTUs were classified against SILVA v132 database [23] and archaeal OTUs against Fungene 8.3 mcrA database [24] with MOTHUR v.1.33.3 [25].

Statistical data analysis

Multivariate analyses were performed in R [30]. Principal coordinate analysis and environmental fitting with packages vegan 2.6-4 [31]. Functional stability was calculated as the coefficient of variation of the main functions of the system. Statistical differences between stability of both reactors were calculated with the asymptotic test for the equality of CVs and the modified signed-likelihood ratio test (SLRT) for equality of CVs using the package cvequality [32]. Sparse Partial Least Squares models, utilized to detect relationships among environmental and biological data, were computed with package MixOmics [33]. ANOSIM and SIMPER analysis were performed in PAST3 [34].

Metagenome annotation and MAGs construction

For shotgun analysis, Trimmomatic software [20] was used for filtering and trimming the raw reads. Metagenome annotation was done in the MG-RAST platform [35] against KEGG [36] and RefSeq databases [37]. Assembly of contigs from reads was done with MEGAHIT [26] and binning with the program MetaBAT [27]. CheckM [28] was used to calculate the completeness and contamination of the bins, meanwhile the taxonomic classification was performed with the program AMPHORA2 [29]. MAGs annotation was made with EnricM software [38] against KEGG [36], pfam [39], and CAZy [40] databases. MAGs Phylogenetic tree was constructed using more than 400 marker genes with PhyloPhlAn [41] and bootstrapping was done with RaxML [42]. Genes significantly enriched in each taxonomic group were determined with t-test (two-sided) implemented in STAMP software [43]. For metabolic modelling and competition/cooperativity indexes MAGs were annotated with Prokka [44] against UniProt database [45]. The KO annotated data were loaded in RevEcoR [46] for metabolic reconstruction and prediction of microbial interactions. Pearson correlations between co-occurrence and metabolic complementarity/competition indexes were computed with function rcorr from hmisc R package [47].

Results

System description and reactor performance

The study involved two anaerobic reactors of different configuration, an UASB and a closed digester tank with a solid recycling system (CDSR), located in the same brewery plant and fed with the same wastewater. Both reactors were sampled 6 times over two years (Supp. Figure 1). Operational and environmental parameters are described in Supp. Figure 1. Despite being industrial wastewater, temperature showed a typical seasonal pattern. Organic loading rate (OLR) and biogas production followed a similar seasonal trend, peaking during summer, when beer production was at its highest. Treatment efficiency remained relatively stable throughout the studied period. Although both reactors received the same wastewater, all measured variables displayed significant differences between reactors (Supp. Table 1), including the functional stability (COD removal efficiency) which was higher for the UASB (t-test, p = 1.9E-4).

Microbial community composition

The bacterial community of 16 samples collected over two years was examined by amplicon sequencing of the V3-V4 hypervariable region of the 16S rDNA, yielding 4,504,974 reads, which after quality filtering, merging, and subsampling totalized 14,250 sequences per sample. Sequences were grouped into 1539 operational taxonomic units (OTUs) at a distance of 0.97 and classified according to Silva v132 [23]. From these sequences, 19% corresponded to Archaea and 81% to Bacteria. Both reactors showed similar bacterial profiles at the phylum level, with the dominance of Proteobacteria (14%-17%), Synergistetes (13%-14%), Bacteroidetes (12%-14%) and Chloroflexi (8% in the CDSR and 18% in the UASB rector), followed by Nitrospirae and Firmicutes (Supp. Figure 2). Within Proteobacteria, 72-78% of the sequences belonged to the class Deltaproteobacteria, orders Desulfuromonadales, Syntrophobacterales, and Desulfovibrionales. Chloroflexi was represented in a range of 80-93% by an uncultured genus of the family Anaerolineaceae, showing 98% of identity with the sequence of the Bacterium JN18_A7_F (DQ168648.1) from a sediment-free PCB-dechlorinating enrichment culture [48].

On the whole, OTUs representing syntrophic bacteria accounted for 18% of total sequences (Synergistales & Syntrophobacterales) plus 13% of Anaerolinaceae, a semi-syntrophic family of bacteria.

The diversity of methanogenic archaea was studied by amplicon sequencing targeting the gene of the alpha subunit of methyl coenzyme M reductase (mcrA), a marker gene. The initial number of reads was 4,479,236. After quality filtering, merging, and subsampling remained 27,224 sequences per sample, which were grouped into 128 OTUs at 85.8% of similarity, a level accepted as equivalent to the 97% for 16S rDNA [22]. At the order level, the community was dominated by Methanobacteriales, with more than 87% of the sequences, followed by Methanomicrobiales (7%). The community composition reflected the prevalence of the hydrogenotrophic pathway of methanogenesis over the acetoclastic metabolism (Supp. Figure 3).

Archaeal communities were highly uneven in both reactors, with only three OTUs accounting for more than 60% of the sequences. For the CDSR, dominant OTUs were classified as Methanobrevibacter arboriphilus (36%), Methanobacterium formicicum (21%), and Methanobacterium beijingense (10%), and for the UASB reactor as Methanobacterium oryzae (27%), Methanobacterium formicicum (23%) and Methanobacterium beijingense (12%).

For comparative purposes, archaeal community composition was also assessed with 16S rDNA amplicons representing 19% of the total sequences. This analysis retrieved 30 OTUs at a distance of 0.03. Dominant archaeal orders according to the 16S rRNA gene were: Methanobacteriales 57%, Methanosarcinales 39%, and Methanomicrobiales 4%. From these, Methanobacterium formicicum was dominant in both reactors with 41-44% of the sequences, followed by Methanosaeta concili with 19-18% of the sequences. This result was highly contrasting with that of mcrA sequencing. Thereafter, subsequent analyses were carried out with mcrA OTUs, which covered a greater diversity among the methanogens of this system.

Bacterial community composition was similar between reactors, implying a strong deterministic effect of wastewater composition on the assembly of the bacterial communities. On the other hand, important differences in the composition of archaeal communities between reactors indicate the influence of other factors on community structure. These results suggest that differences in dominant archaeal members are related to performance of COD-removal function.

Functional diversity in anaerobic communities

Metagenomes of four samples corresponding to the first and last sampling date for each reactor, were analysed by shotgun sequencing. A total of 11.153.303 reads with an average length between 277 and 315 bp remained after quality filtering, from which 82-91% could be classified against the RefSeq database [37]. Microbial community composition was dominated by Bacteria (72% of the sequences) followed by Archaea (28%). According to the MG-RAST automated processing pipeline [35] ca. 83% of sequences matched annotated proteins whilst the rest of the predicted features corresponded to unknown proteins. Of the 35% of unassembled sequences classified against the KEGG Orthology database [36], an average of 60% was related to metabolism. Within this category, 33% were associated with amino acid metabolism, 20% with carbohydrate metabolism, and 16% with energy metabolism. The latter category was dominated by methane metabolism (62%) and oxidative phosphorylation (28%) (Supp. Figure 3). The in-depth analysis of the methane metabolism genes in the systems revealed that 9.3% of the genes were exclusive for hydrogenotrophic methanogenesis, 0.5% were exclusive for acetoclastic methanogenesis, and 0.3% were characteristic of methylotrophic methanogenesis, in addition to 4.0% of genes common to all the methanogenic pathways (Supp. Figure 4). There were no significant differences in functional categories between reactors at any level of classification.

Primary drivers of microbial diversity

The alpha diversity of 16S OTUs assessed by the Shannon-Wiener index or Simpson index showed no significant differences between reactors. On the contrary, the evenness of archaeal OTUs (mcrA) was significantly higher in the UASB reactor (p-value: 0.042; t-test) (Supp. Figure 5). The number of bacterial OTUs was one order of magnitude higher than the number of archaeal taxa, indicating a greater diversity for Bacteria, as clearly illustrated by rarefaction curves (Suppl. Figure 6).

Beta diversity was depicted employing a Principal Component Analysis (PCA) plot based on Bray–Curtis dissimilarity. Bacterial OTUs analysis resulted in samples grouped according to the year of sampling and reactor type, whilst the analysis based on archeal OTUs showed that samples strongly grouped according to the reactor of origin (Figure 1). OLR and week (sampling date) showed a significant correlation with bacterial ordination. Similarly, OLR and VFA partly explained reactor grouping for archaeal communities.

Figure 1. PCA of Archaeal (A) and Bacterial (B) communities representing samples from CDSR (red) and UASB (blue) reactors. Ellipses depict 95% confidence intervals, in violet for samples of 2014 and green for samples of 2015. Arrows indicate significant correlations with environmental factors.

The relative importance of environmental variables (reactor type, operational-physicochemical factors, and sampling date) on the microbial community composition at the OTU level was analysed by multiple regression trees (MRT). The variation in the archaeal community was primarily predicted by reactor configuration (Residual error 0.193). On the contrary, the bacterial composition was partitioned according to a sampling date threshold (week 46, Residual error 0.688). The contribution of reactor type and sampling date to the total variance assessed using the varpart function revealed that reactor type and sampling date made distinct contributions. Reactor configuration accounted for a significant portion of the archaeal community variance (86,5%), while sampling date had a greater impact on the variance in the bacterial community (23,7%).

Analysis of similarities (ANOSIM) was used to confirm the observations made by the MRT analysis. Bacterial composition showed significant differences between years 2014 and 2015 (R: 0.7551; p-value < 0.002). The same analysis applied to the archaeal composition revealed significant differences between reactors (R: 1; p-value < 0.001). A deeper insight into the differences between reactors revealed that abundances of two of the dominant archaeal OTUs account for 49% of the total contribution as was determined by SIMPER analysis. In the CDSR, the community is dominated by an OTU classified as Methanobrevibacter arboriphilus (36% of the mcrA sequences), while the UASB reactor is dominated by a taxon closely related to Methanobacterium oryzae (27% of the mcrA sequences). Regarding bacterial populations, BOTU 4 (Aminiphilus sp.) and BOTU2 (Anaerolineaceae family) prevail during earlier sampling dates, while BOTU7 (Chlorobiales order) and BOTU12 (Desulforomonadales) become dominant after week 56.

Operational parameters as predictors of community structure

For the archaeal community, the environmental predictors of the differences between reactors were VFA and OLR (Figure 1(A)). Specifically, VFA showed greater values for the UASB reactor, while OLR was higher in the CDSR. A rise of inlet COD over time could explain the temporal clustering observed in the bacterial community (Figure 1(B), Supp. Figure, 7). Additionally, OLR appeared as a predictor of differences between bacterial communities in the reactors (Figure 1(B)). The analysis of multiple regressions using the operational parameters as predictors allowed the identification of two major clusters of taxa (Figure 2). The first, starting from the left, is composed of several archaeal OTUs, including two classified as Methanobrevibacter arboriphilus (AOTU1/BOTU10), one closely related to Methanobacterium aarhusense (AOTU11) and two members of the Synergistaceae family (BOTU4 and BOTU 64). This cluster showed a positive correlation with OLR and total alkalinity and correlated negatively with COD-removal efficiency and pH. On the opposite situation, AOTU3 (belonging to Methanobacteriaceae), AOTU12 (classified Methanobrevibacter sp.), BOTU11 (Anaerolineaceae), AOTU4 (classified as Methanospirillum), and BOTU30 (Acidobacteria c5LKS83_ge) were positively correlated with inCOD, pH and higher efficiency in COD removal.

Figure 2. Correlations among community structure and environmental variables calculated by sparse Partial Least Squares (sPLS).

The UASB reactor was starved from week 11 to 51 due to a shortage of brewery production. The effects of starvation were analysed using mixed models, with reactor, period (before, during and after starvation), along with the interaction between reactor and period, considered as fixed factors, and time and interaction time/reactor as random effects. Bacterial diversity (SW) and bacterial evenness (P) showed no significant effects (p > 0.05) from either fixed or random factors. On the contrary, archaeal diversity and evenness showed significant interactions between the reactor and period (p < 0.05 and p < 0.01, respectively), suggesting a potential association with the impact of the starvation period. None of the variables used in the mixed model or their interactions yielded significant effects on the system functioning, assessed as COD removal efficiency (Supp. Figure 8).

Recovery of draft genomes for the dominant members of the prokaryotic community

Draft genomes of the most the dominant OTUs were assembled from the metagenomes of four samples, corresponding to the initial (week 9) and final sampling dates (week 94) from both reactors. The shotgun Illumina sequencing yielded an average of 3.7 10⁶ reads and 8.7 10⁸ bases. After processing, 30 draft genomes were assembled from the metagenomes (MAGs) (Supp. Table 2). Seven MAGs were recovered for the domain Archaea, covering the four most abundant archaeal OTUs (AOTU1, AOTU2, AOTU3 and AOTU5) and the three most abundant 16S rDNA archaeal OTUs (BOTU1, BOTU5 and BOTU13) (Supp. Table 2). For the domain Bacteria, a total of 23 MAGs were assembled, representing the most abundant, but also some rare BOTUs (Supp. Table 2). The OTUs corresponding to six assembled genomes corresponding OTUs displayed significant differences between reactors. BOTU41, classified as Desulfovibrio sp. (MAG I16) and AOTU1, closely related to Methanobrevibacter arboriphilus, were more abundant in the CDSR, whereas AOTU3 (MAG I22), closely related to Methanobacterium oryzae, BOTU11, affiliated to the genus Anaerolinea (MAG R2.19), and BOTU 21 (MAG R2.14) and BOTU 163 (MAG I9), classified as Aminomonas sp., had higher abundances in the UASB reactor.

Functional potential profiles of metagenome-assembled genomes

MAGs were grouped according to their phylogenetic relationships in groups from A to F and annotated against KEGG [36], pfam [39], and CAZY [40] databases (Figure 3). The analysis focused on the differential abundance of annotated genes among phylogenetic groups to elucidate the functional niche occupied by the dominant members of the communities (Supp. Table 3). Group A was composed of four MAGs belonging to Firmicutes, which were significantly enriched in several genes involved in carbohydrate metabolism (Supp. Table 3). Group C exhibited similar abundance in both reactors, and was associated with protein degradation. It consisted of four MAGs, three belonging to the class Deltaproteobacteria, genera Syntrophobacter and Desulfovibrio, and the fourth classified as Thermodesulfovibrio, in the phylum Nitrospirae. Group C was enriched in genes coding for peptidases and branched amino acid transport system, as well as one gene coding for pyruvate-butanoate metabolism. Additionally, it was enriched in several genes implicated in sulphate reduction. Within the Deltaproteobacteria subgroup, the gene coding for alcohol dehydrogenase, as well as two genes related to sulphur metabolism, were detected at high abundance. It is interesting to note that Desulfovibrio genus abundance was higher in the CDSR reactor. Group D included two clusters, one related to Proteobacteria and the other to the Planctomycetaceae family. The MAGs belonging to this group showed an enhanced abundance of genes related to hydrolytic enzymes participating in lipid and protein metabolism. The Planctomycetaceae cluster was enriched in several peptidases and two genes for carbohydrate utilization. The other subgroup includes four MAGs belonging to Proteobacteria, classified as Pseudomonas spp., Serratia and Stenotrophomonas (Gammaproteobacteria), and Bradyrhizobiaceae (Alphaproteobacteria). The Gammaproteobacteria MAGs had a higher level of a serine aminopeptidase gene. Group E, in turn, comprised four MAGs classified as Anaerolinea and Chlorobi, which seemed to be related to carbohydrate and amino acid metabolism. There were only two genes with significantly higher abundance common to this group, one involved in amino acid and propionate degradation, and the other, in carbohydrate hydrolysis. Chlorobi subgroup was enriched in peptidase and malate dehydrogenase genes, meanwhile Anaerolinea subgroup was enriched in glycosyl hydrolases and simple sugar transport system genes. Members of the Anaerolineae subclass showed higher abundances in the CDSR reactor. The group F, include four MAGs classified as Synergistaceae. This group has the potential to exploit sugar, amino acid and fatty acid metabolisms. Group F had significantly higher values of transport-related genes, for C4-dicarboxylates and simple sugars. Additionally, it was enriched in genes related to pyruvate/butanoate metabolism and carbon fixation. Members of the Synergistaceae subgroup were also enriched in several genes involved in glyoxylate and dicarboxylate metabolism, amino acid metabolism and peptide fermentation. Meanwhile, the Aminobacterium cluster showed an increased abundance of a tricarboxylic acid transport gene.

Figure 3. Phylogenetic tree showing relationships among assembled MAGs and the functions significantly enriched in each group. Full circles represent bootstrap support values higher than 98%.

Archaeal MAGs were clustered in group B. All of its members, belonging to Euryarchaeota phylum, displayed genes for hydrogenotrophic methanogenesis. Formate dehydrogenase gene, involved in formate-dependent methanogenesis, was enriched in subgroups B1 (Methanobrevibacter) and B3 (Methanobacterium). Genes required for methanogenesis from methanol and methanogenesis from tri, di, and monomethylamines were present only in Methanobrevibacter and Methanosarcina-related MAGs. Methanobrevibacter and Methanobacterium genera displayed higher abundances in UASB and CDSR reactors, respectively.

According to the functions assigned to each of the differentially abundant MAGs between reactors, CDSR was enriched in genes related to methanol and methylamines methanogenesis, protein degradation, sulphate and alcohol utilization. Conversely, UASB reactor was enriched in genes associated with carbohydrate and lipid degradation, as well as amino acid, fatty acid, and propionate fermentation.

Interspecies interactions as predictors of microbial community structure

With the aim of inferring interactions taking place in the prokaryotic community, a co-occurrence network was constructed based on the Pearson correlation among OTUs with abundance higher than 500 sequences per sample (Figure 4). The network was composed of 6 major clusters, each showing exclusively positive correlations among its members. Each guild (i.e. hydrolytic, acidogenic, acetogenic, and methanogenic) was represented within all the clusters. Negative interactions were seen almost exclusively between clusters C and F. Each of this clusters was characterized by a dominant archaeal member in each reactor, such as Methanobrevibacter arboriphilus in cluster F and Methanobacterium oryzae in cluster C (Figure 4). Most of the negative correlations occurred among several AOTUs classified as M. arboriphilus and M. aarusense in cluster F and M. oryzae in cluster C. There were also negative interactions between M. arboriphilus and two taxa from cluster C classified as Thermanaerovibrio velox. Environmental parameters were included and some of them displayed interesting correlations with OTUs. This is the case of OLR showing positive correlations with methanogenic archaeal OTUs in cluster F and negative correlations with methanogenic archaeal taxa in cluster C. Additionally, pH was positively correlated with BOTU 11 (Leptolinea sp.) and through it with cluster B.

Figure 4. Co-occurrence network showing abundance correlations among OTUs. Red and blue edges correspond to positive and negative relationships. Nodes are coloured according to function; their size is proportional to total abundance. Diamonds depict OTUs significantly more abundant in CDSR and hexagons OTUs significantly more abundant in UASB. OTUs classified as capable of syntrophy have a black outline.

To understand the relative importance of habitat filtering vs species assortment in the community assembly of the studied systems metabolic interaction indexes were calculated by metabolic modelling of the MAGs assembled from the metagenomes. It was found that co-occurrence, measured as Jaccard similarity, positively correlated with the metabolic competition index calculated for the MAGs corresponding to the most abundant taxa (p < 0.01, Pearson correlation). In contrast, co-occurrence correlated negatively with the metabolic complementarity index (p < 0.01, Pearson correlation) suggesting that the metabolic interactions do not define co-occurrence patterns.

Discussion

This study aimed to detect deterministic factors shaping microbial community structure and composition in two full-scale anaerobic reactors with different configurations located at the same brewery plant and treating the same wastewater. Bacterial communities resulted to be similar in composition between both systems. Archaeal communities, on the opposite, were clearly affected by the type of reactor process. Several operational parameters differed between reactors. Among them, OLR and alkalinity were the most important to explain the differences in community structure. To explore the biotic factors acting on the community assembly, a co-occurrence network was constructed using taxa abundances derived from amplicon sequencing. This analysis revealed negative correlations among several OTUs dominant in each system, suggesting that biotic interactions may contribute to differences in archaeal composition. To test this idea, metabolic interaction indices were calculated using metabolic modelling of the MAGs assembled in this study. The comparison between co-occurrence network and interaction indices led us to infer that, in the studied system, species co-occurrence patterns appeared to be shaped more by habitat filtering than by metabolic interactions.

Abiotic factors acting on microbial community assembly

If the feeding composition was the main deterministic factor affecting microbial community structure, we predicted that the prokaryotic communities would be similar in both reactors. It held true for bacterial communities, where differences between reactors followed a similar pattern as temporal dynamics. This is consistent with previous observations that described wastewater composition as a determinant of bacterial community structure, both in activated sludge WWTP and in anaerobic digesters [8,9]. Temporal turnover in microbial community composition are intrinsic to engineered and natural systems and can be attributed to internal or external factors [49,50]. Archaeal communities, on the opposite, resulted clearly affected by the reactor type process, in agreement with the meta-study performed by Leclerc et al., who have reported a major effect of the process type over the nature of the wastewater on the diversity of archaeal sequences [13]. The reactor of origin had an important effect on the structure of the community. the main operational parameters influencing archaeal community structure and composition were related to OLR and alkalinity. Several authors have reported organic loading rate as a determinant of methanogenic community structure both at low F/M [51] and high F/M [52]. In general, as the organic loading rate increases, biogas production initially rises to reach a maximum. However, higher loads beyond this point can cause inhibition and changes in the degradation pathways. For instance, Kong et al. have reported that higher OLR correlates with a switch from acetoclastic to hydrogenotrophic methanogenesis [52]. Alkalinity, in turn, is considered a measure of the system buffering capacity. The UASB reactor, characterized by lower alkalinity, probably maintains an archaeal community more robust to pH changes. Noteworthy, there is a group of OTUs, classified as Methanobacterium sp., Methanospirillum sp., Anaerolineaceae and Dojkabacteria, that correlated positively with the COD removal efficiency while showing a negative correlation with biogas yield, possibly indicating a shift towards biomass production.

Community composition and functional characterization

We have assembled 30 draft genomes corresponding to the majority of dominant species of the system studied. Of these, 16 have a high level of completeness (greater than 80%), low contamination (less than 10%), and large N50, which enable them to be used in future work. They may provide novel genomic information on members of the Bacterial and Archaeal domains, many not yet cultured, as is the case of MAG C22 belonging to the newly discovered Ignavibacteriae SJA-28 lineage. Taxonomic classification of the genomes could be performed down to the species, genus, or family level, according to the information available for their close relatives in the databases. The observed taxa corresponded to those identified by amplicon sequencing and were congruent with the results of other authors in similar systems [53–55].

Based on the results from the metagenomic analysis, the ratio of bacterial/archaeal genes was low, in agreement with the prediction of Beraud-Martinez for anaerobic engineered environments [56]. As might be expected by brewery wastewater composition, the metabolic potential of the anaerobic reactors showed a high proportion of genes related to carbohydrate carbohydrate and amino acid metabolism genes. Energy metabolism genes were predominantly dominated by those involved in methanogenesis. Genes involved in methane production from hydrogen and carbon dioxide were the most abundant, suggesting the dominance of the hydrogenotrophic pathway. This result was consistent with the archaeal members identified through mcrA amplicon sequencing analysis and with the high abundance of syntrophic bacteria assessed by 16S rDNA amplicon sequencing. However, our findings differ from those of earlier 16S rDNA-based studies on the composition of methanogenic archaea in brewery WWTPs, which reported dominance of Methanosaeta sp. [57,58] or Methanosarcinales [59], and implied that the acetoclastic pathway was more important than CO2-reduction pathway [60,61].

The studied bacterial communities were dominated by Proteobacteria, similar to previously reported for the nine full-scale brewery reactors by Werner et al. [62]. The main orders detected within the Proteobacteria class were Syntrophobacterales, Desulfuromonadales, and Desulfovibrionales, all three plausible to maintain close interactions with methanogenic archaea [63–66]. The second phylum in abundance was Synergistetes, dominated by Aminiphilus and Aminivibrio, both genera known to be amino acid fermenters [67,68] as well as other members of the family Sinergistaceae [68,69]. In accordance, five MAGs assembled for this phylum displayed enhanced levels of amino acid catabolic and peptide fermentation genes. Brewery effluent composition is characterized by readily biodegradable organic matter: sugar, soluble starch, ethanol, volatile fatty acids [70] as well as a minor proportion of protein [71]. Regarding the guilds involved in anaerobic digestion, the initial step (hydrolysis) was likely carried out in the study system by Firmicutes, which exhibited higher proportion of genes related to carbohydrate hydrolysis and uptake. Members of Firmicutes, along with Bacteroidetes, were previously identified as responsible for hydrolysis in anaerobic digestion [72]. The phylum Bacteroidetes was represented by Macellibacteroides and SJA-28 (class Ignavibacteria). Although Macellibacteroides fermentans is recognized a carbohydrate fermenter [73], the candidate order SJA-28 has not yet a cultured representative species and there is no definite information about its physiology and metabolism [74,75]. In the present study, an assembled MAG classified as Chlorobi, corresponding to the Ignavibacteriae SJA-28 lineage, contained the full path for propionic acid fermentation, as well as peptidases and amino acid degrading genes. Chloroflexi group and Planctomycetaceae MAGs displayed also high abundance of carbohydrate hydrolysis genes. On other hand, lipases coding genes were more abundant in Plactomycetaceae and Proteobacteria related MAGs. Lipase and esterase lipase activities has been previously reported for several members of Planctomycetaceae [76–78], while lipolytic enzymes have also been described for the remaining four proteobacterial members of the group [79–82]. Fatty acid/dicarboxylic acid anabolic genes were enriched in Synergystaceae MAGs, as has been reported for some members of this family [83,84]. MAGs classified as Pseudomonas and Rhodopseudomonas showed a high abundance of genes for beta-oxidation of fatty acids following previous data for both genera [85,86]. Genes related to protein degradation were present in MAGs belonging to almost all phylogenetic groups. Members of the family Anaerolineaceae are carbohydrate and peptide fermenters [87,88], and have been described as semi-syntrophic organisms [88]. Accordingly, genes for propionate metabolism and amino acid degradation were enriched in MAGs belonging to this lineage. Acetate oxidation genes were enriched in the Clostridiaceae group. Members of this family have been previously reported as candidate SAOB in mesophilic anaerobic reactors [89]. Similarly, Synergistaceae MAGs showed high levels of SAOB marker genes, in agreement with their assignment as SAOB in previous studies [41,90]. Finally, a high proportion of microorganisms belonging to the class Thermodesilfovibrionia, phylum Nitrospirae was detected. The MAG recovered in this study belonging to this class displayed high abundance of genes for sulphate reduction, as well as genes for lactate and propionate syntrophic fermentation, in accordance with the metabolic capacities reported previously for this clade [91]. Hydrogenotrophic methanogenesis genes were present in all methanogenic MAGs, whereas only Methanosaeta and Methanosarcina MAGs showed the capacity for acetoclastic methanogenesis [92,93]. Genes for methanol methanogenesis were detected in both Methanosarcina and Methanobrevibacter MAGs. Members of Methanosarcina are known to use methanol for methanogenesis [92]. However, not all the genomes of Methanobrevibacter associated species display genes for methanol utilization [94]. The use of mono, di, or trimethylamine for methanogenesis was encoded in Methanosarcina MAGs, as was described previously [95]. In the same line, formate dehydrogenase, detected in exclusively hydrogenotrophic MAGs, is indicative of their genetic capacity for formate methanogenesis [95].

Regarding community structure, archaeal communities were less diverse than their bacterial counterparts, regardless of whether they were sequenced using 16S or mcrA primers. A similar observation was reported by other authors [96–98], and could be considered a general rule, taking into account the reduced substrate diversity for methanogenic archaea. A lower functional redundancy for this guild has been explained based on its narrow function [99,100]. Archaeal diversity, as indicated by the Shannon-Wiener index, was similar between reactors. However, the evenness, measured as the Simpson index, was higher for the UASB reactor. According to the intermediate disturbance hypothesis [101], the increased evenness in the UASB reactor could be a consequence of the starvation period, which can be regarded as a moderate disturbance with low frequency. The BACI analysis showed no discernible impact of the starvation on community function (i.e. biogas production and COD removal). This suggests a robust functional buffering capacity despite the relatively low diversity of Archaea.

Distinct archaeal diversity patterns were observed between primers targeting the 16S rDNA and mcrA gene. This lack of agreement was also noted by Wilkings et al., who advocated for the complementarity of both approaches [102]. The main discrepancy between the two amplicons was the importance of the order Methanosarcinales, which was underrepresented in mcrA, likely due to a single mismatch with the commonly used primers designed by Luton [19]. On the other hand, the 16S rDNA gene failed to detect Methanomicrobium and Methanospirillum members and retrieved only one-fourth of the richness captured by mcrA. Based on this, other authors support the preference for mcrA, citing its superior resolution compared with the 16S rRNA gene [103,104].

Biotic factors acting on microbial community assembly

Interspecies interactions could affect both community structure [105] and temporal dynamics [106]. While co-occurrence networks are frequently used to analyse interspecies interactions [107–109], it is important to note that co-occurrence patterns may also arise from niche differences and environmental heterogeneity [110]. The analysis presented in this study showed a common co-occurrence network for both reactors, where the presence of competitive species supports the hypothesis of habitat filtering. Similar results were reported previously for the gut microbiome [111], although contrasting findings have also been reported. For example, Zelezniak found that metabolic dependences modulate co-occurrence in microbial communities [16], The incomplete recovery of less abundant members of the microbial community and the lack of completeness of the assembled MAGs may condition metabolic reconstruction and prediction of microbial interactions. It is interesting to note that negative relationships were not detected within the modules corresponding to the individual reactors, aligning with Blanchet’s et al. remark regarding the challenge of measuring strong negative interactions when one of the species disappears [110]. It seems possible that metabolic interactions influence the ecosystem within each reactor, while differences between them may be explained by habitat filtering. To further test the hypothesis arising from this study and to continue unravelling the importance of microbial interactions in species coexistence and microbial community structure, future experiments using controlled co-culture following a bottom-up approach are essential.

Conclusions

This study supports the hypothesis that wastewater composition exerts a strong deterministic effect on bacterial community composition. On the contrary, the archaeal community exhibited differences in both reactors, with ORL and VFA identified as the relevant abiotic factors explaining these differences. The analysis of the relationship between the co-occurrence network and metabolic interactions between prokaryotic communities of both anaerobic digesters allowed us to conclude that community structuration is mostly determined by habitat filtering. From an operational perspective, correlational data suggested that variations in reactor efficiency were positively correlated with bicarbonate alkalinity and negatively correlated with SS and OLR. Furthermore, the effectiveness of COD removal was positively correlated with the dominance of Methanobacterium oryzae and three syntrophic bacterial taxa.

Availability of data and materials

The datasets supporting the conclusions of this article are available in the NCBI repository, BioProject ID: PRJNA947194 under accession numbers SRR23930591/622 and SRR23949945/8, for amplicon and shotgun raw sequences respectively.

Acknowledgements

We thank staff members of the WWTP. This work was partially funded by FONCyT (PICT 2012-2490, PICT 2018-757 and PICT 2021-CAT-II-0018). L.E. and E.L.M.F are career members of CONICET. T.S.R. was CIN fellow 2015.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by Consejo Nacional de Investigaciones Científicas y Técnicas [L.E. and E.L.M.F. are career members of CONICET]; ANPCyT [PICT 2012-2490, PICT 2018-757, PICT 2021-CATII-001]; CIN [T.S.R. was CIN fellow 2015].