Metabolic functional pro � les of microbial communities in methane production systems treating winery wastewater

Metabolic functional pro les of microbial communities in methane production systems treating winery wastewater Miguel Vital-Jacome (  mvitalj@iingen.unam.mx ) Universidad Nacional Autonoma de Mexico Instituto de Ingenieria https://orcid.org/0000-0003-1263-0168 Julián Carrillo-Reyes Universidad Nacional Autonoma de Mexico Instituto de Ingenieria Germán Buitrón Universidad Nacional Autonoma de Mexico Instituto de Ingenieria


Introduction
Agro-industrial wastes have attracted increasing attention during the last decade because they can be used as raw materials in biotechnological processes, producing energy sources, chemical products, and biomaterials [1,2]. One of these promising raw materials are winery e uents or winery wastewaters, produced during winemaking and characterized by high concentrations of Chemical Oxygen Demand (> 200 g COD L − 1 ), Total Solids (> 50 g TS L − 1 ), and ethanol (> 100 g COD L − 1 ), among other compounds [3]. Anaerobic digestion is an attractive option for treating winery e uents. Methane [4] and other products such as biohydrogen [5] or medium-chain carboxylic acids [6] can also be obtained. Recently, highly concentrated winery e uents were successfully treated to produce methane in single-stage and two-stage (acidogenic and methanogenic) reactors, obtaining better performance in the two-stage con guration and reaching higher productivity in thermophilic than in mesophilic conditions [3,7].
Changes in the performance of anaerobic processes are typically attributed to the process con guration, environment, and operating conditions [8]. However, microbial communities, structure, dynamics, and metabolism deserve further exploration to understand the mechanisms behind the strategies that improve methane production or modify the process stability [9,10]. Molecular tools help to determine the composition and dynamics of anaerobic microbial communities, but understanding the role of the different microbial groups is still one of the biggest challenges [11]. Consequently, several tools were developed to predict functional metabolic pathways and key functional genes based on 16S rDNA data and taxonomic gene information [12]. These tools include PICRUSt and Tax4Fun, which are widely used to study anaerobic digestion processes. For example, Li et al. [13] recently used Tax4Fun to determine the in uence of the substrate-to-inoculum ratio on the metabolic pathways of food waste anaerobic digestion. Bedoya et al. [14] and Gao et al. [15] used Tax4Fun and PICRUSt, respectively, to study the functional pro le of microbial communities in full-scale anaerobic sludge digesters. Zheng et al. [16] used Tax4Fun to study the effect of carbon-to-nitrogen ratios on the metabolic pathways in pig manure and corn stover codigestion. Nevertheless, these studies generally do not relate the process con guration or operating conditions with microbial communities and their metabolism.
This work aimed to compare the taxonomic and functional pro les of four process con gurations of an anaerobic digestion system treating winery e uents. Microbial communities were identi ed using 16S rDNA marker genes, and metabolic functions were predicted using Tax4Fun2, a second-generation tool that includes the prediction of functional gene redundancy. The analysis of these pro les explains, for the rst time, the effect of different process con gurations, such as separated stages (acidogenic-methanogenic), temperature conditions, and changes in operating parameters, such as organic loading rate (OLR) and hydraulic retention time (HRT), over microbial communities and metabolic pathways related to methane production and process stability.

Materials And Methods
Anaerobic digestion system Two con gurations (single and two-stage processes) were compared in the anaerobic digestion system. Mesophilic (35ºC) and thermophilic conditions (55ºC) were applied under four different operational conditions (Table 1). Four conditions were evaluated. Condition 1 consisted of a single-stage methanogenic reactor (M1) operated at mesophilic temperature, a Continuous Stirred Tank Reactor (CSTR) of 7 L (4 L of working volume). Condition 2 consisted of a two-stage process (A2, acidogenic reactor, and M2, methanogenic reactor). Both were operated in mesophilic conditions. The acidogenic reactor was a CSTR of 1 L (0.6 L of working volume), and the methanogenic reactor was a CSTR of 7 L (4 L of working volume). Condition 3 also consisted of a two-stage process (A3-M3), similar to condition 2, but operated at mesophilic (acidogenic) and thermophilic (methanogenic) conditions. Finally, condition 4 consisted of the same con guration and temperatures as condition 3 but with an increment in the organic loading rate (OLR) for both reactors (A4-M4) ( Table 1). The start-up conditions and the acidogenic and methanogenic performances were previously reported [3,7]. The single-stage system (M1) and the acidogenic reactors (A2, A3, and A4) were fed with untreated winery e uents. The e uent of the acidogenic reactors was fed to the methanogenic reactors in the two-stage con gurations (M2, M3, and M4). The acidogenic reactors were subsequently operated from con gurations 2 to 4, i.e., the biomass was not replaced when changing from one con guration to another. The methanogenic reactors of con gurations 1 and 2 were subsequently operated; then, for con guration 3, the biomass was replaced by an inoculum acclimated to thermophilic conditions; con gurations 3 and 4 were subsequently operated. Regarding the comparison among methanogenic reactors, conditions 1 and 2 aimed to compare single and two-stage processes. Con guration 3 was compared with con guration 2 to determine the effect of thermophilic conditions in methanogenesis. Furthermore, con guration 4 was compared with con guration 3 to determine the effect of harsher operating conditions (High-OLR, Low-HRT). ratio (g COD g − 1 COD); f Acidi cation degree (%), i.e., the ratio of the net produced VFA ( nal-initial), expressed as equivalent COD to the initial COD concentration in the reactor feed; g Speci c CH 4 production rate (Nm 3 CH 4 m − 3 d − 1 ); h CH 4 yield (L CH 4 kg COD − 1 ). Data from Vital-Jacome et al. [3] and Vital-Jacome and Buitrón [7].

<Table 1>
Sampling, DNA extraction, and sequencing Samples of 20 mL were taken from each reactor at the end of each operating condition (1 to 4). It is essential to mention that the end of each operating condition was de ned considering the stable operation of each reactor. The criterion for acidogenic reactors was a coe cient of variation of less than 10% in the COD of the acidogenic e uent. A coe cient of variation of less than 10% in methane production was considered for methanogenic reactors. Samples were immediately stored at − 20°C. DNA was extracted from 0.25 g of concentrated biomass of every sample, using the DNeasy PowerMax Soil Kit (Qiagen, Netherlands) and quanti ed by spectrophotometry

Data analysis
For the taxonomic pro les, the sequences were analyzed using the software package "DADA2" v1.6 in the R environment [17]. Brie y, highquality forward and reverse sequences were truncated to -280 and 220-for bacteria and -270 and 200-for bacteria and archaea. Pairing, assembly, and chimera removal were performed on the sequences before taxonomic classi cation using the SILVA release v138 database. Alpha diversity indices, including Chao1, Shannon, and Simpson, were calculated to compare the richness and diversity among samples using the phyloseq package in R. Most dominant genera were taxonomically related to the information of the global MiDAS4 database (Microbial Database for Activated Sludge) and other literature reports, to obtain an indicator on their possible function in the process [18]. The software package Tax4Fun2 was used to predict the functional pro les in the R environment [19]. Brie y, Tax4Fun2 uses the taxonomic pro le to predict the metabolic pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG). The package assigns the ASVs (amplicon sequence variants) to reference sequences in the SILVA database (release v132), and the taxonomic abundance is transformed to the relative abundance of enzymes encoding genes in the KEGG pathways and orthologues. The functional redundancy was calculated for each sample employing the functional redundancy index (FRI) of the Tax4Fun2 package. The FRI refers to the proportion of community members capable of possessing a particular KEGG function and their phylogenetic relationship. For example, a log ratio (FRI A /FRI B ) greater than 0 indicates that a function is more redundant in A. Conversely, a log ratio below 0 indicates that a function is more redundant in B [19].
Results And Discussion

Taxonomic pro le of the acidogenic reactors
Bacterial communities were identi ed in the acidogenic reactors during con gurations 2 to 4 (A2, A3, A4). The bacterial community mainly comprised the phylum Firmicutes in A2, Bacteroidota in A3, and Proteobacteria in A4 (Fig. S1a). The taxonomic pro le at the genus level ( Fig. 1) shows that the predominant genera changed among the con gurations. However, the most abundant genus displayed a relative abundance higher than 70% in all three cases. The genus Sporolactobacillus, predominant in A2, was previously reported in acidogenic reactors devoted to H 2 production, treating highly concentrated substrates, such as molasses wastewater and winery e uents [20,21]. In such systems, this genus was related to the production of lactic acid and volatile fatty acids (VFA). The genus Prevotella, predominant in A3, was reported with relative abundances between 14 to 57.5% in acidogenic reactors treating food waste and sugar re nery wastewater and was related to VFA production, including acetic and butyric acids [22,23]. The genus Acetobacter, predominant in A4, was reported in acidogenic reactors treating tequila vinasses [24]. Acetobacter produces acetic acid from sugars and ethanol-rich environments under anaerobic conditions, so it is commonly found in wine bottles [25], explaining its presence in the treatment of winery e uents. The second most important genus in the acidogenic reactor was Lactobacillus, found in all three con gurations. This genus was identi ed in the fermentation of several substrates, including tequila vinasses [24], molasses [26], cheese whey [27], and it is commonly found in the acidogenic stage of two-stage anaerobic digestion systems treating food waste [28]. In acidogenic reactors, the presence of Lactobacillus is related to the production of lactic acid from carbohydrates before other genera consume the lactic acid to produce VFA [27]. The alpha diversity indexes of community richness (Chao1) and community diversity (Shannon and Simpson) were calculated for the acidogenic reactors (Fig. S1b). The results of these indices indicate a loss of richness and microbial diversity among the acidogenic reactors, which were subsequently operated from conditions 2 to 4.

<Fig. 1>
It was found that the structure of the microbial community changed substantially from con guration 2 to 4, along with a loss of community richness and diversity. However, it should be noted that despite environmental changes, the acidogenic reactors consistently produced acetate as the primary fermentation product in all three con gurations ( Table 1). The main environmental changes in this reactor were related to differences in the OLR and variations in the in uent composition. Most studies have found that variations in OLR and pH drive signi cant changes in the structure of microbial communities of acidogenic reactors. However, despite these community changes, acidogenic reactors can produce a consistent pro le of fermentation products [29]. The loss of richness and diversity in the acidogenic reactors could be bene cial, as Yin et al. recently reported for 42 sludges from acidogenic lab-scale reactors [30]. Those authors suggested that the presence of some critical microorganisms and metabolic functions and not diversity should be considered when optimizing the acidogenic process. As noted in this taxonomic pro le analysis, signi cant changes in the structure of microbial communities do not provide su cient information to understand the response of the process to environmental changes, so it is necessary to analyze the functional pro le of these communities.

Metabolic pathways and functional redundancy in the acidogenic reactors
The functional analysis comprises metabolic pathway categories from the KEGG database distributed in three levels. At the shallowest level (level 3; Fig. S2a), no difference was observed in the acidogenic reactors between con gurations 2, 3, and 4. Also, no difference was observed at level 2 in the metabolism category (Fig. S2b). Meanwhile, Fig. 2 shows the most profound level (level 1) categories for carbohydrate metabolism and energy metabolism, where only slight differences were found. In carbohydrate metabolism, categories such as starch and sucrose metabolism, galactose metabolism, and amino sugar and nucleotide sugar metabolism predicted a slightly higher relative abundance for A2 and A3 than for A4. The main difference in energy metabolism was a higher relative abundance in Oxidative phosphorylation A4 > A3 > A2, probably related to Acetobacter, the most abundant genus in A4, and its ability to grow under aerobic conditions [25]. The present study shows that the abundance of the metabolic pathways predicted for the acidogenic reactors was very similar, despite the signi cant differences in the structure of the microbial community during reactor operation. These results indicate that the microbial community's functional pro le is decoupled from the taxonomic pro le, which is better known as functional redundancy. This concept means that various taxonomically different species share similar or identical functions and can be replaced by one another to maintain the reactor's functionality [31]. Other studies have reported that functional redundancy is an essential feature in acidogenic reactors [32]. Functional redundancy is directly related to the degree of stability and resilience to disturbances in the microbial community [33].
The Venn diagram in Fig. 2c shows that the microbial community within the acidogenic reactors shared up to 5277 metabolic functions during its operation in the three con gurations. Most of these functions are related to the metabolism category. Figure 2c also shows that unique metabolic functions were gradually lost in the acidogenic reactors, going from 348 in A2 to 82 in A3 and only 9 in A4. Since the acidogenic reactors were subsequently operated, this analysis suggests a specialization process of the metabolic functions; most of the unique functions in the community were lost in favor of conserving the shared metabolic functions that allowed the acidogenic process. The functional redundancy of these 5277 shared functions was compared among con gurations. The FRI results showed that more functions had higher functional redundancy in A2 compared with A3, A2 compared with A4, and A3 compared with A4 (Fig. S3). The reduction of the functional redundancy in the acidogenic reactors indicates that functional genes existed in fewer species with a closer phylogenetic relationship at the end of the operation. In general, it was found that in the acidogenic reactors, a loss of diversity was accompanied by metabolic specialization and a reduction in functional redundancy. This nding could indicate the potential destabilization and reduced microbial community resilience. However, the relationship between diversity and functional redundancy in the acidogenic process deserves further study due to insu cient information on this topic and the evidence in recent reports [30].

Taxonomic pro le of the methanogenic reactors
Bacterial and archaeal communities were identi ed in the methanogenic reactors during con gurations 1 to 4 (M1, M2, M3, M4). The most abundant phyla were Firmicutes, Synergistota, Bacteroidota, and Proteobacteria for bacteria, Halobacterota, and Euryarhaeota for archaea (Fig. S4a). The taxonomic pro le, including bacteria and archaea, shows essential differences in the community structure among con gurations (Fig. 3). These differences are discussed in the following subsections, comparing single-stage vs. two-stage methanogenic reactors (M1 vs. M2), mesophilic vs. thermophilic conditions (M2 vs. M3), and changes in operating parameters, from Low-OLR and High-HRT to High-OLR and Low-HRT (M3 vs. M4).

Taxonomic pro le M1 vs. M2
The taxonomic pro les of the methanogenic reactors in con gurations 1 and 2 were compared. It was found that bacteria dominated the microbial community in M1, with functions possibly related to the formation of precursors for producing methane, such as acetate or lactate from Lactobacillus and Citrobacter [34], or H 2 and CO 2 from bacteria possibly belonging to the group of syntrophic acetate-oxidizing microorganisms, such as Syner-01 [35]. In the single-stage reactor (M1), Methanosaeta dominated the archaeal fraction. That archaeon is commonly found in anaerobic digesters performing acetoclastic methanogenesis [18,36]. On the other side, we found that the microbial community in M2 was remarkably more diverse than in M1 (Fig. S4b). The function of the most abundant genera in M2 may be fermentative metabolism, which is directed toward protein fermentation, related to the bacteria genera Proteiniphilum and Thermovirga [18]. A diverse community of archaea carried out methane production in M2. Acetoclastic methanogenesis can be attributed to Methanosaeta, while hydrogenotrophic methanogenesis can be attributed to Methanobacterium, Methanocorpusculum, Methanolinea, and Methanoculleus [18]. According to the diversity indexes (Fig. S4b), M2 had higher microbial richness and diversity than M1. Several studies have compared microbial communities in single and two-stage con gurations. Most studies evidence a higher diversity of bacteria and archaea in the two-stage con gurations [37,38]. Other studies found no differences [39], and some studies reported higher diversity in the single-stage con guration, but they were made under thermophilic conditions [40]. The results of the present study suggest that the twostage con guration promotes microbial diversity in the methanogenic reactor, allowing the growth of microorganisms with a greater diversity of metabolisms, which may explain the superior stability and resilience of the process compared to the single-stage con guration [3].

Taxonomic pro le M2 vs. M3
The differences between the taxonomic pro le of con gurations 2 and 3 (M2 vs. M3) could explain the differences between the process performance in mesophilic and thermophilic conditions. It was found that the archaeal community was more abundant under thermophilic conditions (65%) than under mesophilic conditions (20%). This difference in community structure could explain the higher yields obtained under thermophilic conditions, which were closer to the maximum theoretical yield (Table 1). Another possible explanation could be related to the e ciency of the methane production mechanism between temperature conditions. Under mesophilic conditions, the predominant mechanism in M2 could be a mixture of acetoclastic and hydrogenotrophic methanogenesis carried out by several archaea genera. The microbial community in M3 was composed of Methanothermobacter, a thermophilic hydrogenotrophic methanogen; Syntrophaceticus, classi ed as a syntrophic acetate-oxidizing microorganism; and Acetomicrobium, which is a fermentative producer of acetate, hydrogen, and CO 2 [18]. Therefore, under thermophilic conditions, the predominant mechanism in M3 could be the syntrophic acetate oxidation coupled to hydrogenotrophic methanogenesis, carried out by a single predominant genus of archaea. Regardless of the higher yields, it has been reported that the stability of anaerobic digestion is affected under thermophilic conditions, which is directly related to the microbial community [9]. According to the diversity indexes (Fig. S4b), M2 had higher microbial richness and diversity than M3. Several studies have compared and reported a higher diversity of microbial communities in mesophilic than in thermophilic conditions using substrates such as food waste and sewage sludge [41,42]. This work con rmed the differences between community structures and the higher diversity in mesophilic conditions, which could provide greater resilience and stability to the mesophilic process. Figure 3 shows that changing to harsher operating conditions in the methanogenic reactor (High-OLR, Low-HRT) affected the taxonomic pro le even when the duration was only seven days. For example, the relative abundance of Methanothermobacter increased by 8.1%, and some genera of bacteria disappeared, such as De uviitoga, capable of producing acetate, H 2, and CO 2 as fermentation products. In addition, and according to the diversity indexes, the richness and diversity of the community were signi cantly reduced (Fig. S4b). The reported effect of raising the OLR and reducing the HRT is the loss of communities with particular metabolic functions [43,44]. In the thermophilic methanogenic reactor, bacteria capable of using acetate and other substrates were lost after the harsher conditions. The rst effect of this change in operating conditions was a signi cant increase in methane production. However, after this increase in productivity, it was reported that the process suffered instability and overloading, gradually reducing the process e ciency [7]. Therefore, it can be inferred that applying harsher operating conditions affects the resilience and stability of the process through the loss of microbial communities.

Functional pro le of the methanogenic reactors
According to the KEGG database categories, the metabolic pathways and functions were predicted for the methanogenic reactors in all con gurations (M1, M2, M3, M4). Some of the main differences were related to the Energy metabolism category (level 1; Fig. S5a). The relative abundance of predicted functions related to methane metabolism was twice as high in the con gurations under thermophilic conditions (M3 and M4) compared to those under mesophilic conditions (M1 and M2). This result can be attributed to the higher proportion of archaea in the thermophilic con gurations. Other slight differences in metabolism were observed in the categories of Sulfur metabolism and Carbon xation pathways in prokaryotes. Some other meaningful differences were found in speci c functions of the general categories of Cellular community, Cell motility, Membrane transport, and Signal transduction (level 1; Fig. S5b). The main differences showed a higher abundance in M2 of metabolic pathways related to quorum sensing, bio lm formation, and metabolite transport, among others. These results could indicate that bio lm formation is favored in the two-stage mesophilic methanogenic reactor. However, the level of detail of this rst functional approach does not allow for relating the taxonomic pro le with the functions or mechanisms proposed for each abundant genus of microorganisms.
For this reason, a more detailed analysis was carried out with the results of Tax4Fun2, in which the relative abundance of speci c functions or genes was determined. Figure 4a shows the main identi ed pathways involved in methanogenesis and other routes related to methane production: the methylotrophic methanogenesis pathway (1-4); the hydrogenotrophic methanogenesis pathway (5-11); the acetoclastic methanogenesis pathway (12-16); gene functions common to all methanogenic pathways (17)(18)(19)(20)(21)(22)(23); the Wood-Ljungdahl pathway (24)(25)(26)(27)(28)(29)(30)(31)(32) related to the syntrophic oxidation of acetate [45]. Figure 4b shows the predicted relative abundance of the functions identi ed in Fig. 4a for the methanogenic reactor in the four process con gurations. It was found that the functions related to methane production from methanol and methylamine were only present in M2, probably related to its higher diversity of archaea.

Functional pro le M1 vs. M2
The functional pro le between the single and the two-stage con guration (M1 vs. M2) was similar in all methanogenic pathways. This result indicates that single-and two-stage methane production combines the hydrogenotrophic and acetoclastic pathways under mesophilic conditions in both con gurations. The functions of the Wood-Ljungdahl pathway were also found in M1 and M2. Syntrophic acetate-oxidizing microorganisms use this pathway in reverse (oxidative direction), promoting acetate decomposition into H 2 and CO 2 , supporting the metabolism of hydrogenotrophic methanogens, and competing for acetate with acetoclastic methanogens [45]. The abundance of genes related to this pathway in M1 and M2 could con rm the presence of microorganisms carrying out syntrophic acetate oxidation. As previously suggested, in M1 and M2, this function is probably associated with abundant bacteria genera such as Syner-01. It has been reported that the presence of ammonium, acetate, temperature, retention time, and trace elements in uence syntrophic acetate oxidation [45,46]. The presence of functions of the Wood-Ljungdahl pathway may also be related to the homoacetogenic pathway, which forms acetate from H 2 and CO 2 in the reductive direction [47]. The presence of homoacetogens supports the metabolism of acetoclastic methanogens and competes with hydrogenotrophic methanogens. The abundance of Wood-Ljungdahl pathway functions suggests that some homoacetogens could be present; however, their function could not be related to the most abundant genera in the single-and twostage con gurations under mesophilic conditions. It should be noted that the functional analysis allowed us to support the hypotheses about the methanogenic pathways that yielded the taxonomic information of M1 and M2, but with a higher level of detail.
3.4.2 Functional pro le M2 vs. M3 Figure 4b shows that the functional pro les substantially differ between mesophilic and thermophilic conditions (M1 and M2 vs. M3 and M4). The abundance of genes related to hydrogenotrophic methanogenesis was higher in M3 than in M2, especially the enzyme 4Fe-4S ferredoxin that assimilates H 2 and CO 2 . This result suggests that this pathway is fundamental in producing methane in thermophilic conditions, which is con rmed by the predominance of Methanothermobacter in M3. Interestingly, the enzyme acetyl-CoA synthetase (14) was highly abundant under thermophilic conditions, suggesting that this is the primary mechanism of acetate assimilation, which could be related to the presence of acetoclastic methanogenesis in M3. However, it has been reported that certain enzymes of acetoclastic methanogens are shared by bacteria performing syntrophic acetate oxidation [46]. For this reason, a better explanation for the high predicted abundance of acetyl-CoA synthetase (14) and acetyl-CoA decarbonylase (15,16) is that they are part of the metabolic machinery of syntrophic-acetate oxidizing bacteria. Moreover, the abundance of formate dehydrogenase (30)(31)(32) indicates the presence of this other pathway by forming H 2 and CO 2 from formate.
Another interesting result is the low abundance of other enzymes related to the oxidative Wood-Ljungdahl pathway in thermophilic conditions (26-29). Due to its relative abundance in the taxonomic pro le, the most likely candidate to perform syntrophic acetate oxidation in M3 is Syntrophaceticus. Only one characterized species represent this genus, Syntrophaceticus schinkii, in which acetate oxidation proceeds via the Wood-Ljungdahl pathway [48]. Nevertheless, it has been suggested that syntrophic-acetate oxidizing bacteria can use other metabolic strategies than the Wood-Ljungdahl pathway [49,50], which may be the case of the Syntrophaceticus species identi ed in M3. This hypothesis, the possible alternate pathway followed by syntrophic-acetate oxidizing bacteria, arises from the functional analysis, which would not be possible with taxonomic information alone. However, validating this hypothesis requires the further use of other molecular tools, and further studies are needed to understand the metabolic mechanisms and the ecological role of syntrophicacetate oxidizing microorganisms. Regardless of the mechanism, syntrophic acetate oxidation may replace acetoclastic methanogenesis during the two-stage thermophilic digestion of winery e uents. Previous studies reported that this phenomenon has occurred in other thermophilic methanogenic systems in acetate-rich environments [45,51]. It was recently suggested that a better comprehension of the factors affecting syntrophic acetate oxidation could be used to design strategies to improve the anaerobic digestion processes [47].

Functional pro le M3 vs. M4
In order to nish the analysis of the functional pro les, Fig. 4b shows that the methanogenic reactor's predicted metabolic functions were not signi cantly affected when harsher operating conditions were applied (High-OLR, Low-HRT). This result indicates that despite a loss of diversity and microorganisms genera between M3 and M4, the harsher conditions did not affect the abundance of mechanisms for methane production. However, it is essential to mention that the loss of members of the microbial community could have affected other metabolic functions that were not in the scope of this analysis or even the rate of metabolic functions, which cannot be inferred with taxonomic and functional pro les. Further studies are needed to fully understand the effect of operating parameters on the stability of the process and the metabolism of its communities.

Functional redundancy in the methanogenic reactors
The last stage of the analysis evaluated the functional redundancy in con gurations 1 to 4 for the methanogenic reactors (Fig. S6). A total of 6881 functions are shared among the four con gurations. The overall functional redundancy was higher in the two-stage con guration versus single-stage (M2 > M1); it was higher in mesophilic versus thermophilic conditions (M2 > M3), and it was lower after applying harsher operating conditions (M3 > M4). Fig. S7 compares the functional redundancy of the target metabolic functions related to methane production in the methanogenic reactors. The functional redundancy of these genes was 3 times higher (log ratio of 1.6) in the two-stage con guration compared to the single-stage con guration. Mesophilic conditions provided 2.5 times higher functional redundancy (log ratio of 1.3) than thermophilic conditions. Meanwhile, increased OLR and reduced HRT decreased the functional redundancy up to 1.5 times (log ratio of 0.54) on average. As mentioned above, functional redundancy is directly related to the degree of stability and resilience to disturbances in the process. These results indicate that the two-stage con guration provides higher functional stability and resilience than the single-stage con guration; the microbial community in mesophilic conditions has greater functional resilience than those in thermophilic conditions and applying harsher operating conditions reduces the functional redundancy and stability of the process.

Conclusions
This study showed that the higher stability and resilience of the two-stage con guration were caused by higher microbial diversity and functional redundancy compared to the single-stage. Mesophilic conditions provided higher stability and functional resilience than thermophilic conditions, and applying harsher operating conditions decreased the process stability. Metabolic specialization and a reduction in functional redundancy occurred in the acidogenic reactors, while syntrophic acetate oxidation may play a vital role in the methanogenic reactors under thermophilic conditions. This analysis expands our knowledge of microbial diversity and its relation to microbial functions and process operating conditions, aiming to improve anaerobic process engineering. Figure 1 Taxonomic and functional pro les of the acidogenic reactors. (a) Relative abundance of the predominant genera. "Others" refer to the pool of identi ed genera with a relative abundance lower than 1%. Taxonomic pro le in the methanogenic reactors. Relative abundance of the predominant genera of bacteria and archaea. "Others" refer to the pool of identi ed genera with a relative abundance lower than 1%.