- 1Departmento de Ingenierίa Agroindustrial, Universidad Autónoma Chapingo, Texcoco, Mexico
- 2Department of Biology, West Virginia State University, Institute, WV, United States
- 3Gus R. Douglass Institute, West Virginia State University, Institute, WV, United States
The great diversity of anaerobic digestion (AD) microbiomes indicates high redundancy and flexibility in the assembly of the community. Moreover, AD microbiomes are frequently subjected to disturbances during start-up and operation that require (re)assembly. We tested the reproducibility of secondary succession and AD community assembly mechanisms using a pre-assembled microbiome that was subjected to intense disturbances. Microbiome diversity and functions were followed in replicate mesophilic batch digesters initiated with multiple stressors, including high feed-to-inoculum ratio and many foreign species. Three 10 L batch digesters were derived from a single long-term CSTR digester pre-adapted to poultry litter feedstock and operated in parallel. Physicochemical parameters (methane, acetate, propionate, butyrate, pH, N-NH3, COD) were measured. Metagenome samples were used to assess diversity and functions. Three performance phases were found along the successional gradient: (1) methane inhibition, (2) high methane production, and (3) low methane plateau. The inventory of species (>1600) remained nearly the same, however the relative abundance of species, families, and functions changed during each successional stage. Syntrophic bacteria peaked in abundance during the mid-succession, high methane stage. Succession of overall KEGG functions was highly similar although species and carbohydrate functions diverged during late succession, suggesting diversity of niche partitioning during degradation of recalcitrant organic matter. We estimated the relative contributions of stochastic and deterministic processes and found a shift in the balance during succession. Early succession was not dominated by either dispersal or selection while late succession was dominated by variable selection. In conclusion, methane production recovered following severe (non-lethal) disturbance in a pre-adapted digester microbiome through a reproducible community assembly pathway that shifted toward deterministic, variable selection over time.
1 Introduction
Understanding the drivers and controls of microbial community assembly has become a major theme in microbial ecology (Zhou and Ning, 2017). This is particularly important for engineered microbiomes where the desired functions of the microbiome depend on community structure. Community assembly is the ecological process where microbial species colonize and interact to form local communities through immigration from a regional species pool (Zhou and Ning, 2017). Community assembly can occur in uncolonized environments or following disturbance of an established microbiome. Both contexts represent ecological succession and are found in natural and engineered microbiomes, such as the human gut (Guittar et al., 2019) and waste water treatment plants (de Celis et al., 2022). Originally defined for plant ecology, primary succession occurs in uncolonized environments while secondary succession involves the (re)assembly of communities following substantial disturbances (Prach and Walker, 2011). Conceptually, both types of succession can be applied to industrial microbiomes (Fierer et al., 2010). However, the multitude of bioreactor designs, variations of substrates and products, experimental designs, and different types of stressors, makes studying assembly mechanisms challenging. Because of the importance of controlling assembly in engineered microbiomes, the drivers and controls of succession need to be further understood.
The pathways for microbial community assembly during succession involve both stochastic and deterministic mechanisms (Zhou and Ning, 2017). Stochastic processes are those involving ecological drift (random changes in population abundance) and dispersal. Dispersal can occur at a high rate which can homogenize community composition or low rate which allows for greater influence by drift, and greater community divergence. Deterministic processes, in contrast, are niche-based and emphasize the functional uniqueness of species, including metabolic capacity and interactions (Zhou and Ning, 2017). Selection is the mechanism whereby deterministic processes, coupled to environmental conditions, shape community structure.
The balance between stochastic and deterministic processes can change during community assembly depending on the degree of selection or dispersal, and environmental variability. Four ecological processes reflect shifts in this balance: homogeneous selection, variable selection, homogenizing dispersal, and dispersal limitation (Stegen et al., 2015). Homogeneous selection occurs in environments that are spatially uniform, leading to more similar communities, whereas variable selection occurs in heterogeneous environments that provide greater niche diversity and increased divergence among communities. Homogenizing dispersal occurs when high dispersal rates supersede selection, and interacting communities develop similar compositions. Dispersal limitation occurs when movement between communities is at a low level and selection is also weak. In this case, drift becomes the primary process, leading to divergence of communities. Finally, it may happen that neither dispersal nor selection dominate community assembly; these situations were designated as undominated by Stegen et al., (2015).
Secondary succession occurs following disturbances caused by biotic or abiotic events that cull populations, disrupt community functions, and alter habitats (Shade et al., 2012). Disturbances can be due to stochastic historical events (i.e., contingency) that alter the stochastic/deterministic balance as well as environmental context for community assembly (Zhou et al., 2014). Disturbances can therefore propel succession down multiple pathways that differ in species and functions. Indeed, multiple disturbances may act synergistically resulting in unexpected effects on assembly (Houseman et al., 2008), leading to the development of alternative states where final community compositions are distinct (Fukami and Nakajima, 2011; Fukami, 2015). In engineered microbiomes, the pathways of succession have been shown to be affected by inoculum diversity (De Vrieze et al., 2015; Perrotta et al., 2017), stressors (Regueiro et al., 2015), and stochastic processes (Zhou et al., 2013). This is important because bioreactors are often operated for the products produced by the transient state of the system, rather than final state where metabolic processes change dramatically. Type of disturbance and operational differences, though, make it unclear how reproducible successional trajectories actually are, and difficult to compare diverse systems.
Anaerobic digestion (AD) is a waste treatment process for highly concentrated organic wastes such as livestock manure, producing bioenergy as a byproduct. AD relies on the assembly of a community of bacteria and archaea that requires considerable interspecies cooperation (syntrophy) to convert heterogeneous organic molecules into methane through a multi-trophic level food web. Succession can occur during both digester start-up and operation. Start-up may involve a severe disturbance of a community when inoculum is taken from an established digester. In this case, the substrate to inoculum (S/I) ratio is particularly important for successful start-up and S/I ratios that are too high can lead to acidification and assembly of non-methanogenic communities (Li et al., 2022). Digesters also receive waste feedstock that contains thousands of foreign species that could possibly disrupt community structure, particularly following severe disturbances (Amor et al., 2020; Huber et al., 2020). Nevertheless, long-term operation of digesters with consistent feedstock selects persistent species indicating niche specialization, competitive exclusion, and stability (Werner et al., 2011).
Digesters are typically designed as either flow-through (steady state) or batch systems. Batch digesters represent an interesting case of ecological succession in a closed system. These communities rely on endogenous carbon and energy sources and experience directional change as resources are depleted and metabolites accumulate. However, the primary drivers of succession are not clear. The liquid environment might be expected to mitigate the role of dispersal in community assembly although decomposition of complex particulate organic matter over time could also increase environmental heterogeneity (Allison, 2012; Enke et al., 2018). These factors could shift the balance between stochastic and deterministic processes. In addition, microbial inoculum could be taken from a pre-adapted digester where species sorting and niche partitioning already occurred. Therefore, batch digesters are excellent systems for studying community assembly during secondary succession in bioreactors. Batch digesters are also used in many applications such as Biochemical Methane Potential (BMP) tests (Holliger et al., 2020).
Our objectives were to evaluate the repeatability of community assembly processes in batch mode digesters following disturbances that could produce alternative successional states. If deterministic processes are dominant during succession, it is expected that replicate communities will have highly similar structures due to prior environmental filtering of ecological specialists. However, stochastic processes and historical contingency could lead to alternative successional pathways. Metagenomics was used to investigate whether both species and functional diversity were affected by deterministic or stochastic processes. The start-up of replicate reactors used inoculum from a pre-adapted microbiome that could potentially drive deterministic assembly, but also a very high S/I ratio using animal manure feedstock containing a large reservoir of foreign microbes that could drive stochastic assembly.
2 Materials and methods
2.1 Experimental design
The experiment was designed to test the repeatability of successional pathways and investigate community assembly mechanisms using three replicate batch digesters started with the same pre-adapted microbiome as well as the same feedstock. The inoculum was taken from a 10 L CSTR (continuous stir tank reactor) digester that had been operated long-term with poultry litter. The inoculum was divided equally into three 10 L batch digesters with the same feedstock. The feedstock to inoculum (v/v) ratio was exceptionally high (19:1); far higher than recommended (2-4:1) for digester start-up (Li et al., 2022). The opportunity for stochastic turnover in the community was therefore increased by both feedstock microbial diversity and disturbance due to high dilution. In addition, the use of batch digesters allows a temporal resource gradient to form which facilitates the observation of successional stages and potential shifts in stochastic/deterministic processes.
Three reactor vessels (D1, D2, D3) with 13 L capacity and 10 L working volume were operated for 191 days at 37 ± 2 °C. The feedstock was poultry litter slurry with 3% total solids (TS) from a commercial farm in Tepetlaoxtoc, Mexico. The slurry was made the same day the experiment began by adding water to dry feedstock and filtering to remove large particles. Poultry litter consists of animal excrement, bedding material, and feed residue. The properties of the feedstock were previously reported (Meneses-Reyes et al., 2017). The poultry litter contained TS (66.28%), VS/TS (69.61%), protein (24.70%), ash (20.13%), carbohydrates (19.654%), crude fiber (16.25%), and C:N ratio (8.15) (Meneses-Reyes et al., 2017). The inoculum was 5% (v/v) which came from a CSTR digester adapted to chicken litter for 1283 days. When the inoculum was taken, the parent digester was under steady state conditions (HRT 30 d, 37.6 °C) with pH 7.6, 49.0% methane, and acetate and propionate levels of 1,077 and 682 mg L-1, respectively. Mixing was done once per day by hand rotation.
2.2 Analytical methods
Chemical oxygen demand (COD) was determined using standard methods (APHA standard methods for the examination of water and wastewater, 2005); ammonia nitrogen was determined with HACH Method 10031. Profiles of volatile fatty acids (VFA) were determined as described by (Meneses-Reyes et al., 2018) using gas chromatography (GC) (Claurus 500, Perkin Elmer) equipped with a flame ionization detector and a capillary column Elite-FFAP (length 30 m, diameter 0.32 mm). VFAs were determined using samples of effluent acidified to pH 3 with HCl and centrifuged at 14,500 rpm for 10 min. The supernatant was injected with an auto-sampler using 5 µL volume. Operational conditions were: gas flow 5 mL min-1 at 10.6 psi, 150 °C injection port, 100 °C oven for 8 min with ramp of 160 °C for 8.5 min, 250 °C detector, and retention time of 16.5 min. Helium was the carrier gas. Biogas production was quantified by the displacement method in saline water 10% (v/v) with a flowmeter that uses volume shift and a digital counter (LA8N-BN, Autonics). Methane analysis followed (Guerrero-Toledo et al., 2019) using a GC (Claurus 500, Perkin Elmer). 10 µL biogas samples were injected into the GC with conditions: flow 14 mL min-1 at 14 psi, 100 °C injection port, 70 °C oven, 100 °C detector, retention time 4 min, carrier gas helium. Percentage of methane was obtained using a calibration curve with pure methane. Biochemical Methane Potential (BMP) was calculated following (Meneses-Reyes et al., 2018).
2.3 Metagenomic analysis and statistics
Cell samples were taken three times corresponding to performance phases: beginning (P1), middle (P2), and end-point (P3). To extract cell samples, 50 mL of sludge was pelleted with centrifugation (6700 rpm) and stored at -80 °C. Genomic DNA was extracted using PowerSoil® DNA Kit (Qiagen). Genomic DNA was stored at -80 °C before sequencing. Nine genomic DNA samples were sent to Admera Health (New Jersey) for library construction and sequencing. Library preparation was performed using Nextera XT. Library quality and quantity were assessed with Qubit 2.0 DNA HS Assay (ThermoFisher), Tapestation High Sensitivity D1000 Assay (Agilent Technologies, USA), and QuantStudio ® 5 System (Applied Biosystems). Genomic DNA libraries were sequenced with the Illumina Hiseq 2x150bp format. Sequences were submitted to the National Center for Biotechnology Information Sequence Read Archive under BioProject PRJNA1327443. Sequence quality control used FastQC v.0.11.5. Primers and chimeras were removed using Trim Galore (Martin, 2011) before merging the paired-end sequences using FLASH (Magoč and Salzberg, 2011). All sequences (84.92 Gb) were also deposited in MG-RAST database (Wilke et al., 2015). MG-RAST version 4.0.3 was used for phylogenetic and functional assessments using default parameters (Wilke et al., 2015). Taxon abundance was normalized by dividing by the total number of hits in each metagenome. STAMP version 2.1.3 and R CRAN packages on R Studio and JMP were used for statistical analysis of metagenomic data. Hellinger-transformed data were used. Benjamini-Hochberg FDR correction was used for controlling false discovery rate instead of the familywise error. Scatter plots, hierarchal clustering and Principal Coordinates Analysis (PCoA) were used to evaluate differences in community dynamics and successional trajectories. Physico-chemical parameters (propionate, acetate, butyrate, pH, ammonia nitrogen [N-NH3], chemical oxygen demand [COD], and methane) were evaluated across three operational phases (P1, P2, and P3) using one-way Analysis of Variance (ANOVA). Additionally, Tukey’s Honestly Significant Difference (Tukey HSD) test was used for post-hoc comparisons. Model assumptions were evaluated using residual Q-Q plots and the Shapiro-Wilk test for normality and Levene’s test for homogeneity of variances.
Whether stochastic or deterministic processes were present as community assembly occurred was evaluated using the beta Nearest Taxon Index (βNTI) method (Stegen et al., 2012; Stegen et al., 2013). For this analysis, metagenomic taxonomic profiling was performed with MetaPhlAn 4.0 (Blanco-Míguez et al., 2023) using species-level marker genes. The resulting relative-abundance table was used as the community matrix for ecological distance calculations. A phylogenetic tree corresponding to the detected taxa was used to compute pairwise phylogenetic distances using cophenetic branch lengths with phyloseq function. Taxa absent from the tree or without a matching tip label were removed to ensure consistent taxon sets between the community matrix and the phylogeny. Observed abundance-weighted beta Mean Nearest Taxon Distance (βMNTD) was calculated for each pair of phases using the comdistnt() function in R (Lozupone and Knight, 2005), which quantifies the mean nearest phylogenetic neighbor distance between two communities. This provided a phase × phase matrix describing the observed phylogenetic turnover. Raup-Crick Bray-Curtis (RCbray) analysis was used to partition non-selection processes (drift vs dispersal) using the vegan package (Chase et al., 2011).
3 Results
3.1 Parallel performance of replicate digesters along a successional gradient
Significant differences were observed among operational phases for all measured parameters (methane, propionate, butyrate, acetate, pH, N-NH3, COD) (p < 0.05) indicating distinct process conditions and microbial activity in each successional stage (Table 1, Table 2). Parameters were measured at the same time points as biogas and microbiome sampling. Therefore, three operational phases were defined and associated with stages of community succession. We refer to these phases as: P1, start-up and methane inhibition; P2, maximum methane production; and P3, low methane plateau. Methane concentrations were significantly greater in P2 and P3 than in P1 (p < 0.01), consistent with increased methanogenic activity. During P1, biogas and methane were quite low indicating that dilution with poultry litter feedstock caused inhibition (Figure 1A). During P2, total biogas peaked on day 26 (D1), day 20 (D2), but much later in digester D3 (day 54). As biogas peaked, methane percentage reached 73.2% (D1), 65.5% (D2) and 63.7% (D3). In P3, methane percentage decreased to 31.4% (D1), 25.2% (D2), and 34.8% (D3). COD exhibited a clearly decreasing trend (P1 > P2 > P3, p < 0.001), reflecting progressive organic matter degradation. By P3, COD decreased by 86% in D1, 83% in D2 and 80% in D3 (Table 1). Propionate concentrations were highest in P1, and decreased sharply in P3 relative to P1 and P2 (p < 0.001) indicating enhanced syntrophic oxidation and improved methanogenic conditions (Yidan, 2023). Acetate showed its highest values in P2 (p = 0.018). Butyrate was significantly lower in P3 than in P2 (p = 0.048). pH increased progressively from P1 to P3 (all pairwise comparisons, p < 0.01), coinciding with the consumption of volatile fatty acids and increased ammonia release from protein degradation (Chen et al., 2008). Ammonia (N-NH3) concentrations were significantly higher in P3 compared to P1 (p = 0.003), indicating enhanced deamination activity in later phases. BMP was most similar in D1 and D2. By phase 3, BMP for D3 was 20-29% lower than the other digesters. The progression of metabolic changes was also evaluated with PCA using these performance variables. PCA showed that the phases grouped separately except P2 of D3. PCA component 1 explained 75.2% of the variation while component 2 explained 14.9% (Figure 1). Together, these results indicate distinct metabolic regimes for each operational phase.
Table 2. Summary of one-way ANOVA results for each physico-chemical parameter tested across digester phases (P1, P2, P3).
Figure 1. (A) Biogas and methane production during the three performance phases for each digester (D1, D2, D3). Time points for microbial diversity sampling are indicated by red arrows. (B) PCA analysis of digester performance variables for each successional stage (phase).
3.2 Succession of community structure in replicate digesters
We compared the replicate digester communities to address four questions: Do performance phases correspond to community successional stages? How reproducible are the successional stages between digesters? Do both species and functions display similar trajectories? Are stochastic or deterministic processes dominant during secondary succession in batch digesters? Nine metagenomes were compared representing three performance phases in each digester. Both abundance of taxa and functions were evaluated. Based on classifiable sequences, the relative abundance of Bacteria ranged from 96.1 to 99.01%, and Archaea from 0.5 to 3.38%. Based on MG-RAST, the total taxonomic inventory consisted of 33 phyla, 64 classes, 129 orders, 273 families, 661 genera, and 1622 species. Total identifiable species in each digester was nearly identical: 1602 in D1, 1611 in D2, and 1601 in D3. The majority of species (95%) were shared in every sample.
Overall similarity of digesters was evaluated with multivariate analysis, diversity indices, and correlation analysis. Digester identity showed no detectable effect on community structure with ANOSIM (R = –0.1605, p = 0.7565) and PERMANOVA (F = 0.188, p = 0.8177). The digesters were also compared in terms of Shannon, Simpson, and Fisher’s alpha indices (Table 3). Only Fisher’s alpha had a significant P-value (<0.0005). Fisher’s alpha is sensitive to rare species and high diversity. Using correlation analysis of Bacteria and Archaea separately, all three digesters were quite similar although D3 was somewhat less (Table 4; Supplementary Figure 1). For Bacteria, correlations were: D1/D2, R2 = 0.977; D1/D3 and D2/D3, R2 = 0.943). For Archaea, correlations were: D1/D2, R2 = 0.983; D1/D3, R2 = 0.941; D2/D3, and R2 = 0.945 (Table 4). The overall similarity of total functions between digesters was also very high: D1/D2, R2 = 0.994; D1/D3, R2 = 0.988; and D2/D3, R2 = 0.985. The dominant phyla in P1 were Firmicutes, Proteobacteria, and Actinobacteria; dominance shifted during P2 and P3 to Bacteroidetes, Firmicutes, and Thermotogae (Figure 2A). Abundance of Euryarchaeota varied the most among digesters during P1 (3.3% D1, 1.8% D2, 0.51% D3), but converged by P3 (1.4-1.9%).
Figure 2. Demonstration that digester performance phases represent ecological succession of community structure. (A) Bar graphs showing alpha diversity of the most abundant phyla and families of Bacteria and Archaea during each digester phase. (B) Scatter plots to evaluate changes in the abundance of total species and functions between digester phases (P1, P2, P3). Functions are analyzed in terms of total KEGG functions, energy functions, and carbohydrate functions. For each pair of graphs, the top plot compares P1 and P2, and the bottom plot compares P2 and P3.
Whether performance phases represented community successional stages was considered. Multivariate analysis revealed a strong effect of phases (P1–P3) on microbial community composition. The separation among phases observed in both ANOSIM (R = 0.8025, p = 0.0033) and PERMANOVA (F = 18.136, p = 0.0033) indicates that operational phase exerted a dominant influence on microbial community structure. Scatter plots were used to compare total species in the phases. This showed that P2 diverged considerably from P1 (R2 = 0.300) while P2 and P3 were quite similar (R2 = 0.928), showing again dramatic community turnover (Figure 2B). Family level abundance was also used to evaluate successional gradients because these represent a high degree of shared ecological functions (Goldford et al., 2018; Morrissey et al., 2019). Families that peaked during each phase (p < 0.05) included both high and low abundance groups (Figure 3B). During P1, the dominant families were Lactobacillaceae, Enterobacteriaceae, Corynebacteriaceae and Bacillaceae. The Enterobacteriaceae, in particular, was very abundant (18.8%, D1; 26%, D2; 30.3%, D3) but declined dramatically in P2 (0.6 - 1.6%) and P3 (0.6 - 0.9%). The most abundant species in P1 included, Escherichia coli, Enterobacter sp. 638, Corynebacterium sp., and Salmonella enterica. These species and families are expected in poultry litter feedstock, and are components of chicken guts (Mei et al., 2017; Shang et al., 2018). Many other families also had highest abundance during phase 1 such as Brevibacteriaceae, Burkholderiaceae, and Propionibacteriaceae (Supplementary Figure 2). Dominance during phase 2 (p < 0.05) included families Synergistaceae, Syntrophomonadaceae, Eubacteriaceae, Peptococcaceae, Ruminococcaceae, Thermoanaerobacteraceae, Lachnospiraceae, and others (Figure 3A and Supplementary Figure 2). Family dominance shifted again during phase 3 (p < 0.05) to Bacteroidaceae, Porphyromonadaceae, Thermotogaceae, Prevotellaceae, Flavobacteriaceae, and others (Figure 3A and Supplementary Figure 2).
Figure 3. Box plots showing that particular bacterial families (A), and specific functions (B), change in abundance during the three digester phases (P1, P2, P3) which correspond to three stages of succession. Differences in abundance significant at p <0.05.
The reproducibility of succession between digesters was further evaluated by comparing corresponding phases using scatter plots. For species with relative abundance >0.1%, D1/D2 had highest similarity during P2 (R2 = 0.922); the other comparisons were D1/D3 (R2 = 0.887) and D2/D3, R2 = 0.809 (Supplementary Figure 4). At P3, all three pairings showed decreased similarity, particularly D2/D3 (R2 = 0.573). The less abundant species (0.1-0.02%) were also compared: digesters D1 and D2 were most similar at each phase (Figure 4A). In addition, the level of similarity increased among all three digesters during P2: D1/D2, R2 = 0.873; D1/D3, R2 = 0.807; and D2/D3, R2 = 0.814 (Figure 4A). However, during phase P3, the similarities decreased, and D3 was most divergent: D1/D2, R2 = 0.818; D1/D3, R2 = 0.701; D2/D3, R2 = 0.626 (Figure 2). Beta diversity was measured for species >0.02% using Bray-Curtis dissimilarity and PCoA. This also showed three successional stages (Figure 4C).
Figure 4. Comparison of similarity of each community successional stage (P1, P2, P3) between digesters using scatter plots: (A) species with relative abundance of 0.1 – 0.02% and Bacteroidetes in red, and (B) total KEGG functions. (C) Bray-Curtis PCoA analysis for species with relative abundance ≥0.02%.
Hierarchical cluster analysis was used to evaluate the trajectories of community structure across the successional gradient. This was done by subdividing species into four abundance groups. First, the top 45 most abundant species were considered. The three phases clustered separately and D1/D2 were most similar (Figure 5C). The majority of the abundant species increased during phase 2, including nine species of Bacteroides and five Clostridium species. Second, all species with relative abundances of >0.3% (75 spp.) (Figure 5A) and >0.2% (121 spp.) (Supplementary Figure 5) were compared across all phases. In each case, phases clustered separately and D1/D2 had highest similarity. The most divergent sample was D3 in phase P3. Third, low abundance species (0.1-0.02%) also showed succession: D1 and D2 were more similar to each other than to other phases of the same digester, and, once again, D3 remained distinct in phases 2 and 3 (Figure 5B). Clustering was also used to compare species abundance profiles of several dominant phyla (Firmicutes, Bacteroidetes, Euryarchaeota): D1/D2 showed greatest similarity through all three phases (Supplementary Figure 6). In summary, reproducible successional dynamics were found for both high and low abundance populations, a dramatic shift in species abundances occurred as the methanogenic food web formed, and two digesters were consistently most similar (D1/D2) suggesting two successional pathways might be present.
Figure 5. Hierarchical clustering showing the relative abundance of species: (A) ≥0.3% for all digesters (D1, D2, D3) during all phases (P1, P2, P3); (B) 0.1-0.02% for P2 and P3; and (C) the 45 most abundant species.
To estimate the ecological processes involved in secondary succession we used modeling with βNTI analysis (Stegen et al., 2012). This method determines whether the composition of communities can be associated with random (stochastic) or deterministic processes. If the |βNTI| > 2, deterministic processes are dominant and responsible for the differences between pairwise comparisons of communities. But if the |βNTI| < 2, stochastic processes are dominant (Stegen et al., 2013). In addition, the RCbray metric was used to assess the contributions of drift and dispersal. The three digester phases were compared to each other (Figure 6). The transition from P2 to P3 was classified as deterministic with a βNTI value of 2.60. The other comparisons (P1/P2, P1/P3) had values close to zero indicating a greater contribution of stochastic processes. Following Stegen et al. (2015), the classification of the community assembly processes for each of these transitions is undominated (P1 → P2), variable selection (P2 → P3), and dispersal limitation (P1 → P3) (Figure 6).
Figure 6. Tests of stochastic and deterministic mechanisms during secondary succession in three digesters. Pairwise comparisons were made for the three operational phases (P1, P2, P3) using the digesters as replicates. βNTI values are reported in the boxes. RCbray analysis was used to determine whether dispersal or selection were dominant processes.
3.3 Succession of community functions in replicate digesters
First, we tested whether performance phases represented successional stages. Total functions in P1 compared to P2 showed correlation of R2 = 0.814, while P2 versus P3 was higher (R2 = 0.983) (Figure 2B). This indicates increasing similarity over time. Combined energy functions from all digesters also showed low correlation between P1 and P2 (R2 = 0.630) but much higher correlation for P2 versus P3 (R2 = 0.978) (Figure 2B). The same pattern was found for combined carbohydrate functions (P1/P2, R2 = 0.697; P2/P3, R2 = 0.969) (Figure 2B). Second, to evaluate the reproducibility of successional stages regarding functions, each phase was compared to the corresponding phase in the other digesters. Combined functions showed very high correlations between equivalent phases although D1/D2 was somewhat higher (Figure 4B). Correlations between digesters dropped slightly during P3 (Figure 4B). In addition, two KEGG subsystems were evaluated separately. Comparing total energy functions showed that D1 and D2 were more similar to each other than to D3 during the first two phases (Figure 7A). However, D3 increased in similarity during phase 3. Total carbohydrate functions also had highest correlation between D1/D2 (Figure 7C). In contrast to energy functions, though, carbohydrate function similarity decreased in the final phase (P3), particularly in D3 (Figure 7C). Cluster analysis of carbohydrate and methanogenesis functions also showed clear discrimination between the trajectories of D1/D2 and D3 (Figures 7B, D).
Figure 7. Scatter plots of KEGG energy (A) and carbohydrate (C) functions during the three successional stages (P1, P2, P3). Methane metabolism functions are highlighted in red. Cluster analysis of KEGG methanogenesis functions (B) and carbohydrate functions (D).
Specific examples of energy and carbohydrate metabolism, and cell physiology, were also evaluated separately. Similar to community structure, some KEGG functions peaked (p < 0.05) during each phase. Examples are: phase 1, nitrogen, fatty acid, pyruvate and propanoate metabolisms, phosphotransferase system, and polycyclic aromatic hydrocarbon (PAH) metabolism; phase 2, acetate metabolism (acetate kinase [ackA] and phosphate acetyltransferase [pta] genes) and butanoate metabolism, aminoacyl-tRNA synthetase, ribosomes, and carbon fixation; and phase 3, TCA cycle, glycan degradation, and galactose/fructose/mannose metabolisms (Figure 3B, Supplementary Figure 3).
4 Discussion
We evaluated whether the start-up of batch digesters exhibited repeatable assembly pathways for community structure and function when a pre-adapted methanogenic microbiome was disassembled (diluted) and placed into a similar environment, but with a resource gradient and hundreds to thousands of foreign species. All three replicates exhibited directional change along the gradient which was created by feedbacks between microbiome metabolism and resource quality and quantity. Highly reproducible stages of secondary succession were found that corresponded to three phases of digester performance: (1) methane inhibition (start-up), (2) high methane, and (3) low methane plateau. Each stage was characterized by differences in community structure, including peaks in abundance of particular species and families. We also found successional dynamics for energy and carbohydrate metabolisms and cell physiology. With about 1600 identified species in each digester, the shifts in species and functions during succession were surprisingly similar between digesters. The balance between stochastic and deterministic processes was also found to shift along the successional gradient.
4.1 Stages of secondary succession along a resource gradient
The well-known methanogenic food web contains the steps hydrolysis, acidogenesis, acetogenesis, and methanogenesis. Since methane production requires an integrated food web, both species and functions are expected to reflect successional states. In stage 1 (early succession), the severe dilution of the methanogenic community should disrupt positive interspecific interactions. As expected, methane and total biogas production were temporarily inhibited. However, the restart of methane production did occur fairly quickly indicating recovery of these interactions. To provide H2, CO2 and acetate for methanogenesis, hydrolysis and fermentation of the most labile substrates (sugars, amino acids) are expected to occur first (Pascual-García and Bell, 2020). Supporting this inference, peaks occurred in the phosphotransferase system which transports sugars across membranes (Somavanshi et al., 2016). Propionate levels and propanoate pathway genes were also highest.
Stage 2 (mid-succession) was characterized by peak biogas and methane production and the greatest changes in community structure and function. The resource gradient changed as total organic carbon greatly decreased and food web intermediates (acetate, CO2) increased. This was reflected by a shift in community functions, particularly energy and carbohydrate metabolism. Strong positive interspecific interactions, and high nutrient levels, are expected to underlie maximum methane production (Embree et al., 2015). Evidence for mutualistic interactions involving syntrophic acetate-oxidizing bacteria are peaks in acetate levels, acetate metabolism genes (ackA, pta), and hydrogenotrophic methanogens (Methanomicrobiaceae) (Hattori, 2008; Li et al., 2022). Families containing syntrophic bacteria, Syntrophomonadaceae, Synergistaceae and Thermoanaerobacteraceae, did peak at this time (Zhu et al., 2019). Recolonization of the disturbed environment was further indicated by peaks in ribosome and aminoacyl tRNA biosynthesis genes which represent cell growth (Nemergut et al., 2015). A dramatic change in dominant and minor species also occurred, including the highest convergence among digesters. This was not necessarily expected because of the dilution and high number of foreign species that could potentially invade (Amor et al., 2020; Bittleston et al., 2020). The convergence of minor populations is particularly noteworthy and may represent prior environmental filtering of niche specialists. For example, using replicated freshwater microcosms, Pagaling et al (Pagaling et al., 2014). found greater predictability of assembly when the community was “pre-conditioned” by prior exposure to a particular environment.
Late succession (stage 3) showed the greatest convergence in performance properties. Levels of propionate, butyrate, acetate, COD and pH were most similar among all digesters. Both methane and total biogas substantially decreased, indicating a large decline in readily digestible organic compounds, as well as the lowest levels of propionate and acetate. The original labile substrates are expected to be depleted as decomposition shifts toward recalcitrant organic matter such as cellulose (Rivett et al., 2016). Supporting this hypothesis, carbohydrate functions showed peaks for glycan degradation and metabolism of galactose/mannose which are components of hemicellulose. TCA cycle genes increased while pyruvate metabolism genes decreased indicating shifts from fermentation to respiration. Although high free ammonia can inhibit methanogens, the values here correspond to increased methane production, suggesting that inhibition was minimal or that ammonia-tolerant methanogenic populations were active (Yang et al., 2025). Families peaking included Thermotogaceae, Porphyromonadaceae and Bacteroidaceae which produce enzymes for hydrolysis, as well as polysaccharide catabolism, including cellulose and hemicellulose (Hahnke et al., 2015; Ozbayram et al., 2018; Esercizio et al., 2021). Due to the large reduction of methane and acetate, we predict decreases in strong mutualistic interactions which is consistent with reductions in syntrophic families.
4.2 Stochasticity, determinism, and contingency
We found evidence for a shift in the balance of stochastic and deterministic mechanisms during secondary succession. Using the conceptual system of Stegen et al (Stegen et al., 2015), the transition from stage 1 to 2 can be classified as the undominated process where neither dispersal nor selection were most important. This was likely due to the initial mixing of inoculum and liquid feedstock which created random dispersal in the reactors, followed by community convergence when the methanogenic food web was reassembled. In contrast, the transition from stage 2 to 3 was classified as deterministic with variable selection. We propose that variable selection represented an increase of environmental heterogeneity even though the digester had a liquid environment. This could be due to the slow, irregular breakdown of plant biomass and particulate organic matter in the feedstock (Allison, 2012; Enke et al., 2018).
The increase of deterministic selection during the stage 2–3 transition was associated with divergence of minor populations and carbohydrate functions. We suggest these minor species are specialist degraders of the recalcitrant complex organic matter (Rivett et al., 2016) and utilize carbohydrate-active enzymes (CAZymes) (Berlemont and Martiny, 2016). This is consistent with the expected increase in the heterogeneity of debris biomass during decomposition (Enke et al., 2018). Yu et al., 2024 also found that species with specialized functions can coexist when they share a major resource. These factors may explain why Bacteroidetes, which are specialized polysaccharide degraders (Bertucci et al., 2019), converged in abundance during the initial consumption of labile substrates such as sugars, then diverged during later succession. In addition, the particulate material might have become progressively colonized during the long period of digestion which would further increase environmental variability (Enke et al., 2018). Other possible reasons for the divergence of minor populations could be that the same species have different dynamics in different contexts (Bittleston et al., 2020), or that species interactions (mutualistic versus competitive) changed as nutrients varied (Hoek et al., 2016). Therefore, we propose that the highest degree of determinism was associated with variable selection of minor populations responding to heterogeneous resource quality and spatial complexity.
Deterministic assembly has been indicated in other digester studies as well. For example, Vanwonterghem et al., 2014. found highly parallel changes in community succession in replicate CSTR digesters. However, this study and others (Peces et al., 2018) used quite different experimental conditions, particularly synthetic feedstock. Niche specialization was found by Huber et al., (2020). in a stressed thermophilic digester treating poultry litter when performance recovery coincided with recovery of a specific population rather than replacement by a redundant species. In our study, the parent digester received poultry litter feedstock for more than three years which could facilitate development of mutualistic interactions. The convergence of minor populations also suggests the development of synergistic networks (Nobu et al., 2015; Wu et al., 2016). Previous work with our digesters has shown changes in network organization and modularity during succession which may indicate niche specialists and mutualism (Guerrero-Toledo, 2019). A time-scale dependence for assembly mechanisms was also found. When early succession (P1) was compared to late succession (P3), dispersal limitation was indicated, likely due to the closed nature of the reactors where no new species could immigrate from a regional species pool. This shows that the time scale for sampling successional gradients is important for discerning assembly processes (Dini-Andreote et al., 2015).
Our study also found preliminary evidence for the effects of historical contingency on successional pathways (Fukami, 2015). Although similar successional stages were present in all three digesters, species and functions could be distinguished as two trajectories (D1/D2 and D3) even with the same set of species. For dominant and minor species, as well as functions, the two digesters with the highest initial similarity (D1, D2) showed the most similar successional pathways. In contrast, D3 was less similar and showed delayed production of methane. The difference of D3 at stage 1 was undoubtedly due to a sampling effect. Even though all three digesters were started by equally dividing the inoculum-feedstock mix, D3 showed the highest initial variation in populations. Pascual-Garcia et al (Pascual-García et al., 2025). found that small differences in community composition can lead to different successional trajectories. In our study, the majority of bacteria and methanogens that became dominant in mid and late succession of D1/D2 began in lower abundance in D3 and never reached similar levels. For example, species associated with plant biomass degradation (Clostridium cellulolyticum, C. phytofermentans, C. thermocellum) increased slowly in D3. However, it is important to note that clustering is a descriptive statistic that requires further confirmation, ideally with additional replicate reactors, but this was not feasible given the size and cost of our experiment.
4.3 Limitations and applications
Several limitations and applications of our study can be mentioned. First, an important point is that the disturbed digester microbiomes were already pre-assembled following long-term selection which provided prior opportunity for competitive exclusion and invasion from the feedstock metacommunity. Whether highly reproducible succession would occur with entirely different feedstock is a relevant question. A complementary experiment could use another feedstock with new foreign species and evaluate whether different alternative successional pathways arise. We further speculated that environmental heterogeneity increased in late succession due to progressive decomposition and colonization of organic debris. This hypothesis can be tested by isolating organic particulate matter over time and measuring diversity. Whether alternative transient states are an important feature of AD also deserves further scrutiny. Indeed, one strength of anaerobic digestion is the adaptability of the microbiome toward operational changes suggesting significant breadth in operational community states; nevertheless, normal operation can occasionally slide into dysfunction for no apparent reason. Regarding industrial applications, this study can contribute to defining “microbiome design principles” where controlling bioreactor community assembly is a fundamental objective of the ecological engineering approach (Lawson et al., 2019). Finally, how these findings transfer to steady state digesters remains to be tested. It may be that the undominated mechanism for community assembly recurs with every pulse of fresh feedstock. A deeper understanding of community assembly processes in anaerobic digestion may help to engineer resilience, product quality, and new applications.
Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: NCBI SRA, accession PRJNA1327443.
Author contributions
FG-T: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing. TE-S: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – review & editing. GH-E: Conceptualization, Investigation, Methodology, Resources, Supervision, Validation, Writing – review & editing. DH: Conceptualization, Formal Analysis, Methodology, Resources, Writing – original draft.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. FMGT and TES thank Consejo Nacional de Ciencia y Tecnología (CONACYT) (now SECIHTI), and Universidad Autónoma Chapingo through projects CB-2016-01(288716) and 25017-C-62, respectively, and Leticia Martinez for technical assistance. DHH thanks WVSU Land Grant and USDA NIFA Evans Allen funds (NI241445XXXXG018).
Conflict of interest
The authors declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frmbi.2025.1707779/full#supplementary-material
References
Allison S. D. (2012). A trait-based approach for modelling microbial litter decomposition. Ecol. Lett. 15, 1058–1070. doi: 10.1111/j.1461-0248.2012.01807.x
Amor D. R., Ratzke C., and Gore J. (2020). Transient invaders can induce shifts between alternative sta ble states of microbial communities. Sci. Adv. 6, eaay8676. doi: 10.1126/sciadv.aay8676
Berlemont R. and Martiny A. C. (2016). Glycoside hydrolases across environmental microbial communities. PloS Comput. Biol. 12, 1–16. doi: 10.1371/journal.pcbi.1005300
Bertucci M., Calusinska M., Goux X., Rouland-Lefèvre C., Untereiner B., Ferrer P., et al. (2019). Carbohydrate hydrolytic potential and redundancy of an anaerobic digestion microbiome exposed to acidosis, as uncovered by metagenomics. Appl. Environ. Microbiol. 85, e00895–e00819. doi: 10.1128/AEM.00895-19
Bittleston L. S., Gralka M., Leventhal G. E., Mizrahi I., and Cordero O. X. (2020). Context-dependent dynamics lead to the assembly of functionally distinct microbial communities. Nat. Commun. 11, 1440. doi: 10.1038/s41467-020-15169-0
Blanco-Míguez A., Beghini F., Cumbo F., McIver L. J., Thompson K. N., Zolfo M., et al. (2023). Extending and improving metagenomic taxonomic profiling with uncharacterized species using metaPhlAn 4. Nat. Biotechnol. 41, 1633–1644. doi: 10.1038/s41587-023-01688-w
Chase J. M., Kraft N. J. B., Smith K. G., Vellend M., and Inouye B. D. (2011). Using null models to disentangle variation in community dissimilarity from variation in α-diversity. Ecosphere. 2, art24. doi: 10.1890/ES10-00117.1
Chen Y., Cheng J. J., and Creamer K. S. (2008). Inhibition of anaerobic digestion process: A review. Bioresour. Technol. 99, 4044–4064. doi: 10.1016/j.biortech.2007.01.057
de Celis M., Duque J., Marquina D., Salvadó H., Serrano S., Arregui L., et al. (2022). Niche differentiation drives microbial community assembly and succession in full-scale activated sludge bioreactors. NPJ Biofilms Microbiomes 8, 23. doi: 10.1038/s41522-022-00291-2
De Vrieze J., Gildemyn S., Vilchez-Vargas R., Jáuregui R., Pieper D. H., Verstraete W., et al. (2015). Inoculum selection is crucial to ensure operational stability in anaerobic digestion. Appl. Microbiol. Biotechnol. 99, 189–199. doi: 10.1007/s00253-014-6046-3
Dini-Andreote F., Stegen J. C., Elsas J. D., and van; Salles J. F. (2015). Disentangling mechanisms that mediate the balance between stochastic and deterministic processes in microbial succession. Proc. Natl. Acad. Sci. 112, E1326–E1332. doi: 10.1073/pnas.1414261112
Eaton A. D. and Franson M. A. H. (2005). APHA standard methods for the examination of water and wastewater. 21st edition (Washington, DC: American Public Health Association/American Water Works Association/Water Environment Federation).
Embree M., Liu J. K., Al-Bassam M. M., and Zengler K. (2015). Networks of energetic and metabolic interactions define dynamics in microbial communities. Proc. Natl. Acad. Sci. 112, 15450–15455. doi: 10.1073/pnas.1506034112
Enke T. N., Leventhal G. E., Metzger M., Saavedra J. T., and Cordero O. X. (2018). Microscale ecology regulates particulate organic matter turnover in model marine microbial communities. Nat. Commun. 9, 2743. doi: 10.1038/s41467-018-05159-8
Esercizio N., Lanzilli M., Vastano M., Landi S., Xu Z., Gallo C., et al. (2021). Fermentation of biodegradable organic waste by the family thermotogaceae. Resources. 2, art24. doi: 10.3390/resources10040034
Fierer N., Nemergut D., Knight R., and Craine J. M. (2010). Changes through Time: Integrating Microorganisms into the Study of Succession. Res. Microbiol. 161, 635–642. doi: 10.1016/j.resmic.2010.06.002
Fukami T. (2015). Historical contingency in community assembly: integrating niches, species pools, and priority effects. Annu. Rev. Ecol. Evol. Syst. 46, 1–23. doi: 10.1146/annurev-ecolsys-110411-160340
Fukami T. and Nakajima M. (2011). Community assembly: alternative st able states or alternative transient states? Ecol. Lett. 14, 973–984. doi: 10.1111/j.1461-0248.2011.01663.x
Goldford J. E., Lu N., Bajić D., Estrela S., Tikhonov M., Sanchez-Gorostiaga A., et al. (2018). Emergent simplicity in microbial community assembly. Science 361, 469–474. doi: 10.1126/science.aat1168
Guerrero-Toledo F. (2019). Dinamica Microbiana de Biodigestores Enriquecidos con Propionato. PhD Thesis (Texcoco, Mexico: Universidad Autonoma Chapingo).
Guerrero-Toledo F., de M., Espinosa-Solares T., Balagurusamy N., Guerra-Ramirez D., Huber D. H., et al. (2019). Potencial bioquimico de metano de pollinaza adicionada con propionato en condiciones mesofilicas. Sci. Agropecu. 10, 307–311. doi: 10.17268/sci.agropecu.2019.02.17
Guittar J., Shade A., and Litchman E. (2019). Trait-based community assembly and succession of the infant gut microbiome. Nat. Commun. 10, 512. doi: 10.1038/s41467-019-08377-w
Hahnke S., Maus I., Wibberg D., Tomazetto G., Pühler A., Klocke M., et al. (2015). Complete genome sequence of the novel porphyromonadaceae bacterium strain ING2-E5B isolated from a mesophilic lab-scale biogas reactor. J. Biotechnol. 193, 34–36. doi: 10.1016/j.jbiotec.2014.11.010
Hattori S. (2008). Syntrophic acetate-oxidizing microbes in methanogenic environments. Microbes Environ. 23, 118–127. doi: 10.1264/jsme2.23.118
Hoek T. A., Axelrod K., Biancalani T., Yurtsev E. A., Liu J., and Gore J. (2016). Resource availability modulates the cooperative and competitive nature of a microbial cross-feeding mutualism. PloS Biol. 14, 1–17. doi: 10.1371/journal.pbio.1002540
Holliger C., Astals S., de Laclos H. F., Hafner S. D., Koch K., and Weinrich S. (2020). Towards a standardization of biomethane potential tests: A commentary. Water Sci. Technol. 83, 247–250. doi: 10.2166/wst.2020.569
Houseman G. R., Mittelbach G. G., Reynolds H. L., and Gross K. L. (2008). Perturbations alter community convergence, divergence, and formation of multiple community states. Ecology 89, 2172–2180. doi: 10.1890/07-1228.1
Huber D. H., Ramirez-Garcia A., Chavarria-Palma J. E., Espinosa-Solares T., Lhilhi Noundou V., Montenegro-Garcia N. A., et al. (2020). Stress induced by crude glycerol in a thermophilic digester: microbial community divergence and resilience, but slow process recovery. Appl. Microbiol. Biotechnol. 104, 10769–10781. doi: 10.1007/s00253-020-10965-6
Lawson C. E., Harcombe W. R., Hatzenpichler R., Lindemann S. R., Löffler F. E., O’Malley M. A., et al. (2019). Common principles and best practices for engineering microbiomes. Nat. Rev. Microbiol. 17, 725–741. doi: 10.1038/s41579-019-0255-9
Li Y., Chen Z., Peng Y., Huang W., Liu J., Mironov V., et al. (2022). Deeper Insights into the Effects of Substrate to Inoculum Ratio Selection on the Relationship of Kinetic Parameters, Microbial Communities, and Key Metabolic Pathways during the Anaerobic Digestion of Food Waste. Water Res. 217, 118440. doi: 10.1016/j.watres.2022.118440
Li C., Hao L., Lü F., Duan H., Zhang H., and He P. (2022). Syntrophic acetate-oxidizing microbial consortia enriched from full-scale mesophilic food waste anaerobic digesters showing high biodiversity and functional redundancy. mSystems 7, e00339–e00322. doi: 10.1128/msystems.00339-22
Lozupone C. and Knight R. (2005). UniFrac: A new phylogenetic method for comparing microbial communities. Appl. Environ. Microbiol. 71, 8228–8235. doi: 10.1128/AEM.71.12.8228-8235.2005
Magoč T. and Salzberg S. L. (2011). FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27, 2957–2963. doi: 10.1093/bioinformatics/btr507
Martin M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12. doi: 10.14806/ej.17.1.200
Mei R., Nobu M. K., Narihiro T., Kuroda K., Sierra J. M., Wu Z., et al. (2017). Operation-Driven Heterogeneity and Overlooked Feed-Associated Populations in Global Anaerobic Digester Microbiome. Water Res. 124, 77–84. doi: 10.1016/j.watres.2017.07.050
Meneses-Reyes J. C., Hernández-Eugenio G., Huber D. H., Balagurusamy N., and Espinosa-Solares T. (2017). Biochemical methane potential of oil-extracted microalgae and glycerol in co-digestion with chicken litter. Bioresour. Technol. 224, 373–379. doi: 10.1016/j.biortech.2016.11.012
Meneses-Reyes J. C., Hernández-Eugenio G., Huber D. H., Balagurusamy N., and Espinosa-Solares T. (2018). Oil-extracted chlorella vulgaris biomass and glycerol bioconversion to methane via continuous anaerobic co-digestion with chicken litter. Renew. Energy 128, 223–229. doi: 10.1016/j.renene.2018.05.053
Morrissey E. M., Mau R. L., Hayer M., Liu X.-J. A., Schwartz E., Dijkstra P., et al. (2019). Evolutionary history constrains microbial traits across environmental variation. Nat. Ecol. Evol. 3, 1064–1069. doi: 10.1038/s41559-019-0918-y
Nemergut D. R., Knelman J. E., Ferrenberg S., Bilinski T., Melbourne B., Jiang L., et al. (2015). Decreases in Average Bacterial Community rRNA Operon Copy Number during Succession. ISME J. 10, 1147–1156. doi: 10.1038/ismej.2015.191
Nobu M. K., Narihiro T., Rinke C., Kamagata Y., Tringe S. G., Woyke T., et al. (2015). Microbial dark matter ecogenomics reveals complex synergistic networks in a methanogenic bioreactor. ISME J. 9, 1710–1722. doi: 10.1038/ismej.2014.256
Ozbayram E. G., Ince O., Ince B., Harms H., and Kleinsteuber S. (2018). Comparison of rumen and manure microbiomes and implications for the inoculation of anaerobic digesters. Microorganisms 6, 15. doi: 10.3390/microorganisms6010015
Pagaling E., Strathdee F., Spears B. M., Cates M. E., Allen R. J., and Free A. (2014). Community history affects the predictability of microbial ecosystem development. ISME J. 8, 19–30. doi: 10.1038/ismej.2013.150
Pascual-García A. and Bell T. (2020). Community-level signatures of ecological succession in natural bacterial communities. Nat. Commun. 11, 2386. doi: 10.1038/s41467-020-16011-3
Pascual-García A., Rivett D. W., Jones M. L., and Bell T. (2025). Replicating community dynamics reveals how initial composition shapes the functional outcomes of bacterial communities. Nat. Commun. 16, 3002. doi: 10.1038/s41467-025-57591-2
Peces M., Astals S., Jensen P. D., and Clarke W. P. (2018). Deterministic mechanisms define the long-term anaerobic digestion microbiome and its functionality regardless of the initial microbial community. Water Res. 141, 366–376. doi: 10.1016/j.watres.2018.05.028
Perrotta A. R., Kumaraswamy R., Bastidas-Oyanedel J. R., Alm E. J., and Rodríguez J. (2017). Inoculum composition determines microbial community and function in an anaerobic sequential batch reactor. PloS One 12, 1–14. doi: 10.1371/journal.pone.0171369
Prach K. and Walker L. R. (2011). Four opportunities for studies of ecological succession. Trends Ecol. Evol. 26, 119–123. doi: 10.1016/j.tree.2010.12.007
Regueiro L., Spirito C. M., Usack J. G., Hospodsky D., Werner J. J., and Angenent L. T. (2015). Comparing the inhibitory thresholds of dairy manure co-digesters after prolonged acclimation periods: part 2 – correlations between microbiomes and environment. Water Res. 87, 458–466. doi: 10.1016/j.watres.2015.05.046
Rivett D. W., Scheuerl T., Culbert C. T., Mombrikotb S. B., Johnstone E., Barraclough T. G., et al. (2016). Resource-dependent attenuation of species interactions during bacterial succession. ISME J. 10, 2259–2268. doi: 10.1038/ismej.2016.11
Shade A., Peter H., Allison S. D., Baho D., Berga M., Buergmann H., et al. (2012). Fundamentals of microbial community resistance and resilience. Front. Microbiol. 3 doi: 10.3389/fmicb.2012.00417
Shang Y., Kumar S., Oakley B., and Kim W. K. (2018). Chicken gut microbiota: importance and detection technology. Front. Vet. Sci. 5-2018. doi: 10.3389/fvets.2018.00254
Somavanshi R., Ghosh B., and Sourjik V. (2016). Sugar influx sensing by the phosphotransferase system of escherichia coli. PloS Biol. 14, 1–19. doi: 10.1371/journal.pbio.2000074
Stegen J. C., Lin X., Fredrickson J. K., Chen X., Kennedy D. W., Murray C. J., et al. (2013). Quantifying community assembly processes and identifying features that impose them. ISME J. 7, 2069–2079. doi: 10.1038/ismej.2013.93
Stegen J. C., Lin X., Fredrickson J. K., and Konopka A. E. (2015). Estimating and mapping ecological processes influencing microbial community assembly. Front. Microbiol. 6. doi: 10.3389/fmicb.2015.00370
Stegen J. C., Lin X., Konopka A. E., and Fredrickson J. K. (2012). Stochastic and deterministic assembly processes in subsurface microbial communities. ISME J. 6, 1653–1664. doi: 10.1038/ismej.2012.22
Vanwonterghem I., Jensen P. D., Dennis P. G., Hugenholtz P., Rabaey K., and Tyson G. W. (2014). Deterministic processes guide long-term synchronised population dynamics in replicate anaerobic digesters. ISME J. 8, 2015–2028. doi: 10.1038/ismej.2014.50
Werner J. J., Knights D., Garcia M. L., Scalfone N. B., Smith S., Yarasheski K., et al. (2011). Bacterial community structures are unique and resilient in full-scale bioenergy systems. Proc. Natl. Acad. Sci. 108, 4158–4163. doi: 10.1073/pnas.1015676108
Wilke A., Bischof J., Gerlach W., Glass E., Harrison T., Keegan K. P., et al. (2015). The MG-RAST metagenomics database and portal in 2015. Nucleic Acids Res. 44, D590–D594. doi: 10.1093/nar/gkv1322
Wu L., Yang Y., Chen S., Zhao M., Zhu Z., Yang S., et al. (2016). Long-term successional dynamics of microbial association networks in anaerobic digestion processes. Water Res. 104, 1–10. doi: 10.1016/j.watres.2016.07.072
Yang J., Zhang J., Du X., Gao T., Cheng Z., Fu W., et al. (2025). Ammonia inhibition in anaerobic digestion of organic waste: A review. Int. J. Environ. Sci. Technol. 22, 3927–3942. doi: 10.1007/s13762-024-06029-1
Yidan J. (2023). Lu yahai syntrophic propionate oxidation: one of the rate-limiting steps of organic matter decomposition in anoxic environments. Appl. Environ. Microbiol. 89, e00384–e00323. doi: 10.1128/aem.00384-23
Yu X. A., McLean C., Hehemann J.-H., Angeles-Albores D., Wu F., Muszyński A., et al. (2024). Low-Level Resource Partitioning Supports Coexistence among Functionally Redundant Bacteria during Successional Dynamics. ISME J. 18, wrad013. doi: 10.1093/ismejo/wrad013
Zhou J., Deng Y., Zhang P., Xue K., Liang Y., Nostrand J. D. V., et al. (2014). Stochasticity, succession, and environmental perturbations in a fluidic ecosystem. Proc. Natl. Acad. Sci. 111, E836–E845. doi: 10.1073/pnas.1324044111
Zhou J., Liu W., Deng Y., Jiang Y.-H., Xue K., He Z., et al. (2013). Stochastic assembly leads to alternative communities with distinct functions in a bioreactor microbial community. mBio 4, 10.1128/mbio.00584–12. doi: 10.1128/mbio.00584-12
Zhou J. and Ning D. (2017). Stochastic community assembly: does it matter in microbial ecology? Microbiol. Mol. Biol. Rev. 81, 10.1128/mmbr.00002–17. doi: 10.1128/mmbr.00002-17
Keywords: community, anaerobic digestion, deterministic, microbiome, secondary succession, stochastic
Citation: Guerrero-Toledo FdM, Espinosa-Solares T, Hernández-Eugenio G and Huber DH (2026) Community assembly following disturbance in batch anaerobic digesters displays highly reproducible secondary succession and a shifting stochastic-deterministic balance. Front. Microbiomes 4:1707779. doi: 10.3389/frmbi.2025.1707779
Received: 17 September 2025; Accepted: 18 December 2025; Revised: 10 December 2025;
Published: 27 January 2026.
Edited by:
Muhammad Ali, Trinity College Dublin, IrelandReviewed by:
Hari Ananda Rao, Curtin University, AustraliaYutong Liu, The Pennsylvania State University (PSU), United States
Copyright © 2026 Guerrero-Toledo, Espinosa-Solares, Hernández-Eugenio and Huber. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: David H. Huber, aHViZXJkaEB3dnN0YXRldS5lZHU=
Flor de María Guerrero-Toledo1,2,3