Unveiling the deterministic dynamics of microbial meta-metabolism: a multi-omics investigation of anaerobic biodegradation

Abstract

Background

Microbial anaerobic metabolism is a key driver of biogeochemical cycles, influencing ecosystem function and health of both natural and engineered environments. However, the temporal dynamics of the intricate interactions between microorganisms and the organic metabolites are still poorly understood. Leveraging metagenomic and metabolomic approaches, we unveiled the principles governing microbial metabolism during a 96-day anaerobic bioreactor experiment.

Results

During the turnover and assembly of metabolites, homogeneous selection was predominant, peaking at 84.05% on day 12. Consistent dynamic coordination between microbes and metabolites was observed regarding their composition and assembly processes. Our findings suggested that microbes drove deterministic metabolite turnover, leading to consistent molecular conversions across parallel reactors. Moreover, due to the more favorable thermodynamics of N-containing organic biotransformations, microbes preferentially carried out sequential degradations from N-containing to S-containing compounds. Similarly, the metabolic strategy of C18 lipid-like molecules could switch from synthesis to degradation due to nutrient exhaustion and thermodynamical disadvantage. This indicated that community biotransformation thermodynamics emerged as a key regulator of both catabolic and synthetic metabolisms, shaping metabolic strategy shifts at the community level. Furthermore, the co-occurrence network of microbes-metabolites was structured around microbial metabolic functions centered on methanogenesis, with CH₄ as a network hub, connecting with 62.15% of total nodes as 1st and 2nd neighbors. Microbes aggregate molecules with different molecular traits and are modularized depending on their metabolic abilities. They established increasingly positive relationships with high-molecular-weight molecules, facilitating resource acquisition and energy utilization. This metabolic complementarity and substance exchange further underscored the cooperative nature of microbial interactions.

Conclusions

All results revealed three key rules governing microbial anaerobic degradation. These rules indicate that microbes adapt to environmental conditions according to their community-level metabolic trade-offs and synergistic metabolic functions, further driving the deterministic dynamics of molecular composition. This research offers valuable insights for enhancing the prediction and regulation of microbial activities and carbon flow in anaerobic environments.

Video Abstract

Introduction

Microbial growth and activity drive the global carbon cycle, and an essential part of carbon flow is mediated by microbial anaerobic metabolism [1,2,3]. Approximately, 70% of global methane (CH₄) emissions are created via microbial anaerobic metabolism in typical anaerobic habitats, including the deep layers of terrestrial soils and aquatic sediments, guts and rumens, and anaerobic facilities [4]. Anaerobic microorganisms also contribute to carbon storage via microbial carbon sequestration [3, 5], which leads to the transformation of organic carbon from labile to recalcitrate molecules in terrestrial and marine carbon pools [5,6,7]. Moreover, the anaerobic metabolic activities of animal digestive flora are relevant to host health [8]. Understanding microbial anaerobic metabolism promotes the prediction of and response to global climate change, the regulation of carbon flow during microbial engineering, and the maintenance of vital systems’ health [7, 9, 10]. However, the core rules underlying anaerobic metabolic patterns and dynamics remain largely understudied within natural habitats and engineered systems [11,12,13].

Microorganisms are capable of complex and flexible metabolic pathways [14, 15]. It has been emphasized that microorganisms exhibit different metabolic behaviors to adapt to habitats and resist disturbance [16, 17]. Meanwhile, microorganisms also have metabolic discrepancies in their preferences for substances and forms of metabolic utilization during various activities such as growth, multiplication, and other trophic activities [2, 3]. This has led to controversial inferences and opposing experimental evidence in judging the effects of environmental and/or biotic factors on microbial metabolisms. For example, there are debates regarding the positive and negative correlations of microbial carbon use efficiency with microbial carbon storage [7]. Furthermore, the observed changes in microbial metabolites represent the collective metabolic performances of individual microorganisms under specific environmental conditions [16, 17]. Therefore, we introduce the concept of meta-metabolism to represent the collective metabolic interactions and dynamics of various microorganisms within a given environment. An exploratory study by Danczak et al. [18] has suggested that the composition of metabolites, although not performing active behaviors like microorganisms, can be considered similar to microbial communities, undergoing turnover and assembly, and fulfilling ecological functions. The advancement of multi-omics detection technologies provides vast amounts of data about the complex compositions of microbial communities and metabolites [19, 20]. However, a new challenge arises in linking the complex compositions of microbial communities and metabolites [21, 22]. To address this challenge, it is crucial to characterize metabolic patterns and continuously monitor the response of microbial meta-metabolism to the surrounding conditions [16, 23, 24].

The anaerobic digestion (AD) system can serve as a model for studying the dynamics and rules of microbial meta-metabolism. It operates as a relatively self-contained and closed system characterized by rapid rates of biological metabolism [25, 26], allowing for the observation of the dynamic pattern within the system. As a successful bioengineering technique, AD facilitates the conversion of macromolecular organic matter into intermediates such as volatile fatty acids and biogases [4]. In this resource recovery process, the regulation of carbon flows towards target products is quite essential to promote the efficiency [12, 27]. During the AD process, organic molecules are mainly generated through the ex vivo modification and in vivo turnover of microorganisms as metabolites and released into the liquid environment [3]. The microbial and metabolic compositions within the AD system are incredibly complex [12, 28] and synchronously influenced by microbial cooperation, metabolic thermodynamics, and variable environmental factors [29,30,31]. These co-factors contribute to the diverse metabolic pathways, trade-offs, and strategies within the system, offering a tractable model system for studying the microbial metabolic mechanisms.

Here, we combined metagenomic and metabolomic data to discover the dynamic mechanisms of microbial anaerobic meta-metabolism during the simulated AD experiment. Based on identifying dissolved organic matter (DOM) metabolites, the inference of biotransformation, and the linkages with microbial activities, we focused on the overall and specific dynamic metabolic patterns under anaerobic conditions. The hypotheses of this study were as follows: (i) assembly processes can elucidate the consecutive turnover of DOM metabolites, akin to microbial community dynamics, with harmonious interactions between metabolite assemblages and microbial communities; (ii) microbial communities exhibit metabolic strategy trade-offs throughout the process, which can be understood using general principles derived from meta-metabolism evidences; and (iii) the selection of metabolites by microorganisms depends on microbial ecological functions and molecular traits, influenced by specific environmental conditions.

Materials and methods

Experiment and sampling

Food waste (FW) was selected as the substrate of our bioreactor due to its high volume as a major source of municipal organic waste and a vital application scenario for AD [32]. To simulate the anaerobic digestion of FW and ensure maximum degradation of organic matter [33, 34], we set up six parallel bioreactors that were run continuously for 96 days under anaerobic conditions. The simulated FW used in the study (see Table S1 for details) was modified from a previous study [35]. The FW suspension was prepared by mixing simulated FW and tap water in a 1:3 wet-weight ratio. Each reactor (total volume 2.87 L) was ultimately filled with 2.5 L of the initial starting mixture composed of seed anaerobic sludge and FW suspension in a 1:4 wet-weight ratio. After set up, these reactors were controlled within a temperature range of 35–38 °C [26] and a stirring speed of 100 rpm. Samples were collected on days 0, 3, 6, 9, 12, 24, 48, and 96 after the experiment was set up. Physicochemical properties were measured according to our previous study [25], and more details are provided in the supplementary information (SI).

Metabolomic preprocessing

The untargeted metabolomics proceeded using a Fourier-transform ion cyclotron resonance mass spectrometer (FT-ICR-MS) with a 15.0 T superconducting magnet and a negative electrospray ion source (Bruker solariX, Bruker, USA). Before detection, the DOM metabolites in the filtered liquor (0.45-µm filter) were enriched from the aqueous phase to the organic phase using solid-phase extraction with a Bond Elut-PPL straight barrel (1 g, 6 mL Agilent, USA) according to a previous study [36]. The mass spectrometric detection was performed in negative ion mode, and the signal was scanned and accumulated three hundred times to improve quality [37]. After the quality control process (see SI for details), 11,086 peaks were assigned to distinct molecular formulas. Diverse characteristics of each formula were calculated using R software v4.0.4, including O/C ratio, H/C ratio, N/C ratio, P/C ratio, S/C ratio, Kendrick defect based on CH₂ (kdefect_CH2), double bond equivalent (DBE), DBE minus oxygen (DBE-O), aromaticity index (AI), modified AI (AI_mod), DBE minus AI (DBE-AI), the nominal oxidation state of carbon (NOSC), standard Gibbs free energy of carbon oxidation (GFE), and carbon use efficiency (Y_met) [18, 38]. In addition, analyzed formulas were classified based on their elemental composition compared to natural organic matter [38,39,40]. The molecules were then classified into seven categories: lignin-like, lipid-like, unsaturated-hydrocarbon-like, condensed-aromatic-like, protein-like, carbohydrate-like, and tannin-like (Table S2).

Metagenomic preprocessing

Following the manufacturer’s instructions, microbial DNA was extracted from the time series samples (filter cake obtained by filtering 4 mL of sludge, n = 48) using the PowerSoil™ DNA Isolation® Kit (MO BIO Laboratories, USA). The V4 region of the 16S rRNA gene of prokaryotes was amplified using the universal primer pair 515F and 806R [41]. The procedures of amplification and purification were consistent with our previous study [42]. Sequencing was finally performed on an Illumina NovaSeq platform (Magigene Biotech, Guangzhou, China). The sequenced data were analyzed on our amplicon sequencing data analysis pipeline (https://dmap.denglab.org.cn) [43]. Subsequently, 8319 zero-radius operational taxa units (ZOTUs) were obtained for further analysis. In addition, DNA extracted from each sample was subjected to shotgun sequencing on an Illumina NovaSeq platform (Magigene Biotech, Guangzhou, China), while the mixed DNA from six replicates at each time was subjected to nanopore metagenome sequencing (BENAGEN, Wuhan, China). Metagenome-assembled genomes (MAGs) were then obtained from hybrid metagenomic assembly (see SI for details). The predicted genes of different MAGs were reannotated via eggNOG-mapper online [44], and pathway reconstruction was conducted via KEGG Mapper [45]. Furthermore, DRAM [46, 47] and dbCAN3 [48] were utilized to provide extra gene annotation related to short-chain fatty acid (SCFA) metabolisms and carbohydrate-active enzymes (CAZYmes). More details are provided in the SI.

Diversity and relative abundance analysis

Alpha diversity was calculated and displayed based on observed richness. The dissimilarity between compositions of different samples was computed using the Bray–Curtis distance, and the results were visually represented through principal coordinates analysis (PCoA). The relative abundances of various categories were measured according to the relative peak intensity and the elementary composition of the detected formulas. When analyzing and visualizing the relative abundance of organic matter with varying carbon content, we ranked the average accumulated abundance of each category over time and displayed only the top five in the final results.

Temporal conversion characteristics analysis

The conversion characteristics of metabolites throughout the entire process were observed to infer the metabolic strategies employed. The overall characteristics of the metabolites were calculated based on the weight of the corresponding characteristics of all formulas in each sample [40]. The variations of these general characteristics during the process were used to analyze the trends of metabolite conversion. Then, potential biotransformations [49] were determined based on pairwise mass distance [12, 50] mapping to a professional database [18] (https://github.com/danczakre/Meta-Metabolome_Ecology). The specific biotransformations of different time points were then identified using the linear discriminant analysis effect size (LEfSe) method [51]. If the relative abundance of special biotransformations at a specific time point was significantly higher than at all other time points (p < 0.05), we considered these transformations time-specific. Furthermore, community biotransformation thermodynamics were calculated to elucidate trade-offs in microbial metabolic strategies. The biotransformation thermodynamics of pairwise molecules were determined by calculating the difference in their GFE [49]. A smaller GFE difference indicates that the transformation between paired molecules occurs more easily. To determine the direction of a transformation, which is necessary for calculating GFE difference, we used the increase and decrease in molecular weight between paired molecules. An increase in molecular weight corresponds to synthetic metabolism, while a decrease corresponds to catabolic processes. The mean values of the biotransformation thermodynamics for all pairwise molecules with the same characteristics were defined as the community biotransformation thermodynamics for a specific type of molecular transformation.

Assembly process analysis

Organic matter turnover could arise as the collective result of the gain, loss, and transformation history of individual molecules [18]. This process is similar to the turnover of microbial communities, allowing the extended application of assembly processes to DOM. After that, deterministic processes, influenced by environmental and/or biotic factors, drive the composition to follow a convergent or divergent pattern of assembly, as “homogeneous selection” and “heterogeneous selection,” respectively [52]. Alternatively, the DOM component could be structured according to stochastic processes such as drift and dispersal in microbial community ecology [18, 52]. The inferred community assembly mechanisms, via phylogenetic bin-based null model analysis (iCAMP) [53], were utilized to quantify the process of assembling metabolite composition. Under this framework, metabolites were assigned to a tree based on their similarities of molecular traits, analogous to the phylogenetic tree of microorganisms. The distances between metabolites on the tree were determined by their combined Euclidean distances across 16 molecular trait dimensions (Table S3) [52]. Based on the distribution of metabolites on the tree, and following the default procedures of the software, closely related metabolites are clustered into bins, ultimately forming 57 bins [52, 53]. The molecular assembly process was determined for each bin, and the overall process was weighted according to the total relative abundance of metabolites in these bins.

To assess the relationship between the turnovers of the metabolite composition and the microbial community, Mantel tests based on Spearman correlation were utilized to calculate the associations between these two components, specifically in terms of taxonomic composition variation and ecological niche structure. To accomplish this, we used Jaccard and Bray–Curtis distances to represent differences in taxonomic composition. Additionally, we considered microbial phylogeny and molecular similarity as corresponding concepts and calculated the β-nearest-taxon index (βNTI) to illustrate differences in ecological niches [18, 53, 54].

To evaluate the relative importance of co-turnover of the microbial community and individual environmental variables in metabolite composition, multiple regression on matrices (MRM) analysis was utilized. The procedure involved filtering environmental explanatory variables with higher correlation (Spearman’s ρ² < 0.5) [55]. The matrices of metabolite composition and microbial community were calculated using Bray–Curtis distance, while matrices of environmental explanatory variables were structured based on Euclidean distance. Before the regression, all input distance matrices were standardized, with 0 as the mean and 1 as the standard deviation. Therefore, the coefficients of the regression could be used to approximate the degree of the influence of different explanatory variables. MRM was performed using the R package ecodist v2.0.9 with 100,000 permutation tests. We also repeated the same procedure independently to explain microbial community turnover.

Co-occurrence network analysis

To further explore the relationships among individuals in the system, a bipartite network was constructed between microbes and substances, including DOM metabolites and environmental variables. First, individuals which appeared on half of the sampling days and in at least five of the six replicates for individual sampling days were filtered and retained. The selected individuals were then standardized to a range of 0 to 1, and the standardized data was used to compute local similarity (LS) scores to assess the correlation of the time-series data [56, 57], allowing for a delay of 1-time point. The random matrix theory (RMT) method was employed to select a threshold based on the LS score matrix, and the threshold was set at 0.819 [58]. After module identification, we calculated the relative within-module degree and the participation coefficient of the individuals in the network to determine their role in the network [58]. The network analysis process was performed on the iNAP platform [59] (https://inap.denglab.org.cn).

To evaluate the variation of molecular traits among different network modules, a multi-response permutation procedure (MRPP) analysis was performed based on molecular similarity distance. To analyze the relationship between the network structure and the metabolite’s molecular weight, the local positive correlation ratio was calculated step by step in a molecular weight range of 120, starting from the minimum molecular weight in the network and with a step of 10. Furthermore, to explore the metabolic function of key species in the bipartite network, these species were related to MAGs via local alignment search based on the 16S rRNA gene and annotation comparison.

Results

Discovering intrinsic traits of the overall DOM metabolites composition

Meta-metabolites offer a vast functional space for microbial metabolism, serving as a shared material pool wherein thousands or more metabolites are produced, exchanged, and assimilated by microorganisms [60, 61]. The initial step should be comprehending the composition of this DOM metabolites pool. Using a high-resolution mass spectrum, a total of 11,086 molecular formulas of DOM were identified throughout the 48 samples collected at eight distinct time points (days 0, 3, 6, 9, 12, 24, 48, and 96) in all AD bioreactors. As shown in Fig. S1, many formulas displayed lignin-like (4290) and lipid-like (1,572) compositions. When focusing on the elemental composition of these formulas, it was shown that the majority of analyzed matters contained the element oxygen (O) (10,566/11,086) (Fig. 1A). Consequently, the presence or absence of O was not considered separately in the elemental composition analysis. In descending order, the prevalent categories within this AD system included CH(O) (3149), CH(O)N (2663), CH(O)NS (1749), and CH(O)S (1533) (Fig. S1). In contrast, only a tiny proportion of the identified DOMs (1992/11,086) contained the element P, categorized as CH(O)NP, CH(O)P, and CH(O)PS. The molecules/formulas with 18 C atoms were the most numerous, while the number of molecules gradually decreased with variations from 18 (Fig. 1A). A similar trend was observed when considering H and O elements, where the highest quantity of molecules was reached with 6 H and 22 O (Fig. 1A). In addition, the molecules did not display a significant bias towards unsaturation, with most DOM molecules containing reducing C (Fig. 1B). Therefore, the overall composition of the metabolites in this AD system for food wastes was determined to be biodegradable.

Fig. 1

General characteristics of molecules and the overall variations throughout the AD process. A The number distribution of molecules with different numbers of atoms for C, H, O, N, S, and P. B The degree of unsaturation and redox state of all formulas. (DBE-O)/C and NOSC were used to evaluate these two indexes, respectively. C Variation of alpha diversity over time. Observed richness was used to quantify alpha diversity, pairwise difference was calculated via Wilcoxon test, and p-value was adjusted using Benjamini–Hochberg method. To indicate statistical significance, different letters were assigned to the results. Groups that do not share any letters are significantly different (p < 0.05), while groups that share at least one letter are not significantly different. D Stadardized richness of different categories over time. The results were stadardized into the range 0 to 1 before visualization. E Average relative abundance of different categories across the whole process. F Average relative abundance of the most abundant five components based on number of C atoms. The standard error of the data are represented by error bars in D, E, F. Abbreviations in D and E CB, carbohydrate-like; CA, condensed-aromatics-like; LG, lignin-like; LP, lipid-like; PT, protein-like; TN, tannins-like; UH, unsaturated-hydrocarbons-like; UC, unclassified

The composition of metabolites underwent a continuous and rhythmic succession. Specifically, the observed richness of metabolites increased during the initial stage (days 0–3), followed by fluctuations and decreases throughout the rest of the process (Fig. 1C). This indicates that most metabolites were generated through biological processes accompanied by hydrolysis. However, despite increased richness during days 0–3, we observed a trend towards a composition of simpler metabolites over time (days 3–96). Meanwhile, the continuous turnover of metabolites was visualized using principal coordinates analysis based on Bray–Curtis dissimilarity (Fig. S2). All bioreactors exhibited a consistent pattern in the succession of metabolite composition, wherein the system gradually converted towards a simpler molecular makeup over time. Furthermore, while most of the categories identified in this study followed the trend of overall observed richness, lipid-like substances showed a distinct variation (Fig. 1D). Unlike most substances that only increased at the initial phase (days 0–3), lipid-like substances increased in richness during days 0–9 (Fig. 1D), accompanied by a steady increase in relative abundance until day 12 (Fig. 1E) when methane generation began (Fig. S3). Therefore, the production of lipid-like components in this system exhibited a prolonged phase, during which most substances were hydrolyzed and then other categories decreased. We further analyzed the five most abundant classifications at each time point by carbon chain length. Our results showed that the C18 component was the most-produced throughout the experimental period (Fig. 1F and S4). Moreover, we observed an apparent turnover in the length of the most abundant carbon-chain classes. For example, the relative abundance of the C18 component increased during the first 12 days of the experiment, followed by a subsequent decrease from day 12 to day 96. Meanwhile, other carbon-chain classes that were initially dominant gradually gave way to different classes throughout the experiment. These results indicated the metabolic conversion obeyed definite dynamic rules during the AD process.

Inferring microbial meta-metabolic strategies from temporal characteristics of the metabolites conversion

In addition to revealing variations in metabolite composition, we sought to explore molecular conversion patterns throughout the dynamic process, shedding light on microbial meta-metabolic strategies. Weighted average compositional values were calculated to evaluate the general trend of molecular conversion in the process. We found the overall molecular weight was reduced over the process (Table 1), indicating that anaerobic microorganisms reduced the total amount of organic matter and converted larger molecules into smaller ones. In particular, we observed a sharp decrease in C and O during the early stages (days 0–12) of the process, while C and H primarily decreased during the later stages (days 9–96) (Table 1 and Fig. S5). We also found the degree of unsaturation, redox state, and carbon use efficiency showed a decreasing trend followed by an increasing trend around day 12 (Table 1). Lower values of unsaturation and redox state suggested a higher lability of molecules, and the lower values of carbon use efficiency suggested a higher thermodynamic efficiency of metabolic reactions involved in biomass accumulation [62]. Therefore, our findings indicated that early microbial meta-metabolism might render the entire organic matter composition more susceptible to biodegradation before methane generation. Meanwhile, we observed that the molecular weight of the lipid-like component increased until day 9 (Table S4). This increase in molecular weight of the lipid-like component directly demonstrated their microbial synthesis and elongation during the early stage of the experiment.

Table 1 Weighted average compositional values of DOM metabolites throughout the AD process

After that, we conducted a detailed examination of specific transient biotransformations to discover microbial metabolic actions and trade-offs. Pairwise mass distances were analyzed for each sample, and the results were mapped to 866 common biochemical transformations (Figs. 2A and S6). We further measured their variance over the process to confirm the time specificity of these potential biotransformations. Given that the relative abundance of the majority of potential biotransformations (616/866) at a specific time point was significantly higher than at all other time points (LEfSe, Wilcoxon test’s p < 0.05), we considered these transformations to be time-specific (Fig. 2B). The metabolic transformations performed by the microbial community can be seen as manifestations of their metabolic strategies. Thereby, the differential metabolic actions suggest that the microbial community adopts distinct metabolic strategies throughout the experiment. When considering both the time specificity of the transformations and the elemental composition involved, we found that the transformations of nitrogen (N) and sulfur (S) elements exhibited distinct temporal patterns in this anaerobic biodegradation system. Specifically, N-containing biotransformations were concentrated in the early stages of the experiment (days 0, 3 and 6), while S-containing biotransformations were more prominent in the later stages (after day 12). Subsequent attempts were made to calculate biotransformation thermodynamics at the community level to explore the underlying metabolic strategies. Considering the overall process of molecular biodegradation within the system, we calculated the thermodynamic changes associated with molecular transformations of by subtracting the GFE of smaller molecules from that of larger molecules. This calculation is performed in the direction of decreasing molecular weight. The average value of element-related transformations for each sample represents the thermodynamic difficulty of these transformations within that sample. The analysis revealed that N-containing biotransformations were thermodynamically easier than S-containing biotransformations during microbial anaerobic decomposition (Fig. 2C), suggesting a prioritization of available resources by microorganisms. Meanwhile, the biodegradation of phosphorus (P) containing molecules did not show time specificity during the AD process, which might be typical for microbial metabolic activities with low thermodynamic barriers (Fig. 2C). The switch from synthesis to degradation of C18 lipid-like molecules provided additional evidence that thermodynamic conditions might determine the selection and execution of metabolic strategies by microorganisms. By considering the synthesis of lipid-like substances, the GFE change was then calculated in the direction of increasing molecular weight. We found the thermodynamic barriers to the synthesis of C18 lipid-like substances from small molecules showed significant differences between the synthesis (days 0–9) and catabolism phases (days 12–96) of C18 lipid-like substances (Fig. 2D). These findings on meta-metabolism indicated that microorganisms execute metabolic strategies in response to specific physical conditions such as system thermodynamics and available resources.

Fig. 2

Transient characteristics and thermodynamics of metabolite conversion based on biotransformation analysis. A Molecular biotransformation with time specificity on day 96. Molecular biotransformations were represented using mass distances. The difference analysis was performed using LEfSe, and the time specificity on day 96 was exhibited as an example. The compositions of molecular biotransformations were annotated using different colors. B Time specificity of molecular biotransformation and their elemental compositions. The compositions and time specificities were annotated using different colors. C The thermodynamics of N-containing, S-containing, and P-containing biotransformations during the anaerobic decomposition process. The mean values of the energy change for pairwise molecules as degradation reactions were calculated for each sample. D The thermodynamics of the synthesis of C18 lipid-like substances from small molecules. The mean values of the energy change for pairwise molecules as synthesis reactions were calculated for each sample. The pairwise difference was calculated by the Wilcoxon test in C and D, and the label “***” indicates the p-value < 0.001

Extending community ecology to understand microbial meta-matabolism

Community ecology has provided us with new perspectives on the assembly of metabolite composition [18]. Referring to microbial community assembly [54], homogeneous and heterogeneous selection can imply that metabolite composition becomes convergent or divergent under the influence of biotic and abiotic factors [18, 52]. Homogenizing dispersal and dispersal limitation can indicate whether molecules easily diffuse and migrate within the environment, thereby forming molecular composition [18, 52]. The remaining parts (drift and others) can suggest that molecules undergo random transformations, unobserved processes, or unexpected events that result in the observed molecular composition [52, 54]. We found that homogeneous selection dominated the assembly of metabolites in the system (Fig. 3A). Nevertheless, a pattern of increasing and then decreasing proportion of homogeneous selection was observed (with a peak at 84.05% on day 12), while the proportion of drift and others continued to increase, achieving a level similar to that of homogeneous selection at day 96 (36.53% vs. 36.76%) (Fig. 3A). As the overall ecological process was weighted according to the processes of 57 bins (Fig. S7) of similar molecules, we also found the largest number of individual bins were subject to homogeneous selection during days 0–24, with a range of 21 to 31 bins (Fig. 3A). However, on days 48 and 96, more bins were influenced by drift and others, with a total of 23 and 27 bins respectively, out of which 18 and 16 were driven by homogeneous selection (Fig. 3A). These results suggested the influence of the environment and biotic processes on the organic matter might have a dominant role in determining the content of the large to small molecular transformation. In contrast, an eventual increase in stochastic processes might represent a weakening of the directed conversion activity of microorganisms.On the basis of homogeneous selection being dominant during the molecular turnover process, we next explored the relative attribution and mechanism of microbiological processes. Firstly, the Mantel test affirmed the co-turnover in the compositions of molecules and microbes (p < 0.05, Table S5). To further assess the extent of the drive on molecular turnover from environmental variables and biotic procession, we conducted multiple regression on matrices analysis to explore linear relationships between the response distance matrix and multivariate explanatory distance matrices. When fitting the turnover of the molecular composition, we found that microbial community structure had the most significant influence (coef. = 0.327, p = 0.00001), indicating that microbiological processes played a crucial role in shaping the composition of molecules (Fig. 3B). Similarly, we observed a significant impact of the molecular composition on the succession of the microbial community (coef. = 0.388, p = 0.00001). Additionally, pH significantly influenced microbial community structure (coef. = 0.419, p = 0.00001), suggesting the microbial community was structured (Fig. S8) by both the physical environment (Table S6) and the composition of organic carbon sources.

Fig. 3

Assembly mechanism of DOM metabolites and the correlation with microbial community. A Assembly processes of metabolite component across the anaerobic digestion process. The similarity of molecular properties was calculated based on 16 molecular traits, and then the proportion of different assembly processes was calculated based on the molecular similarity relationships and relative abundance of the metabolite component using iCAMP software. All molecules were clustered into 57 bins during this process. The variance of relative abundance and properties of the assembly process for each bin are shown individually. Each bin is labeled with the molecular formula of the most dominant matter. B The extent of influence of metabolite components and microbes from each other and environmental factors. The extents of influence were evaluated based on the coefficient from the multiple regression on matrices (MRM). The collinearities among environmental factors were obtained based on the Spearman method, and the correlations between individual factors and the turnover of metabolite components and microbes were obtained using the Mantel test based on distance matrices. C The relationship of the similarity of metabolite component and phylogenetic similarity of the microbial community. The correlation was calculated based on the Mantel test using the Spearman method

Following the confirmation of co-turnover in composition between the microbial community and molecular composition, we analyzed the correlations between their ecological functions. Phylogenetically closely related microorganisms tended to share similar ecological niches, forming the foundation for studying microbial ecological processes based on phylogenetic trees [63]. Similarly, molecular chemical traits such as size, availability, and energy provision influenced their specific behaviors in biochemical processes, including synthesis, degradation, and diffusion [18, 21]. A reasonable inference was that molecules with similar molecular traits in multiple dimensions (see the 16 molecular traits in Table S3) exhibited analogous ecological functions. These inferences provided the basis for studying the ecological functional relationships between microorganisms and molecules using microbial phylogenetic trees and molecular trait trees. Mantel test based on Spearman correlation showed a significant correlation between the βNTI values of the two components (p = 0.001) (Fig. 3C), indicating their interdependence in ecological functions. As βNTI was also used to infer the assembly process, the correlation between the βNTI values of the two components also suggested the coordination of the assembly process between microbial community and DOM composition. This implies that the two components underwent similar processes in most cases. Our study found that when one component of microbes and metabolites was subject to selection (|βNTI|> 2), the proportions of the other component structured under the same process were 78.43% (βNTI > 2, χ² = 514.152, p = 0.000) and 43.89% (βNTI <  − 2, χ² = 265.606, p = 0.000) (Table S7). For the stochastic process (|βNTI|< 2), the probability of the consistency of the other component was 24.48% (χ² = 7.567, p = 0.006) (Table S7). These results suggest that the assembly mechanisms of the two components influenced each other due to the correlation of their ecological functions. In particular, each component is a crucial factor of the other during the selection process, coordinating the co-turnover between the microbial community and molecular composition.

Linking microbial functions with metabolites traits

After confirming the mutualistic relationship between microbes and DOM metabolites in the turnover and assembly processes, we further structured the co-occurrence network to infer the underlying relationship between the two components. We treated microorganisms as one group and organic molecules and environmental variables as a separate group to obtain a robust time series network between the two groups. We found that the microbes-metabolites network of the system was assembled around CH₄ generation, with CH₄ identified as a network hub (Figs. 4A and S9). The network was tightly connected, with 62.15% of the total nodes being the first and second neighbors of CH₄. Expanding to the fourth-order neighbors allows connection to almost all nodes (99.49%) (Fig. 4A). Meanwhile, all connector hubs in this network belonged to organic matters (Figs. 4A and S9). Microorganisms provided the second network hub (ZOTU 756, family Ruminococcaceae) and the majority of module hubs in the network (Figs. 4A and S9). Different metabolic activities might be modularly distributed, as evidenced by the detection of significant differences in molecular traits in various modules (Fig. 4B). In addition, we found that larger molecules were more likely to have positive correlations with microbes (Fig. 4C). These monotonically positive correlations with microbes indicated a more unstable position for larger molecules in the network. In general, the system was structured around CH₄ generation, with different microbial metabolic activities aggregating microorganisms and substances into different modules to form the network structure and drive the conversion of molecules into smaller matters.

Fig. 4

Co-occurrence network analysis between microbes and analyzed substances throughout the AD process. A Bipartite network structured around methane generation. The analyzed substances in the network included DOM metabolites and other detected environmental substances. Each node was annotated with its role in the network, its category, and the neighborhood relationship with CH₄. B Difference analysis of organic molecules in different modules based on molecular trait similarity. The distance matrix of molecules was calculated based on 16 traits of the molecular formula. The pairwise difference was quantified using MRPP, and those marked with a star indicate that the difference was significant (p < 0.05). C Local positive correlation ratio in the network varied with molecular weight. Each point represents the ratio of positive correlation with the microorganisms of a series of molecules in a 120-molecular-mass range. D Functional annotation of the genomes of key microorganisms identified in the bipartite network

On the basis of the key species identified in the network, we further matched these key species to their MAGs to associate their key roles with their metabolic functions. Through the comparison of 16S rRNA gene sequences and annotations, seven ZOTUs were successfully matched to four MAGs (Table S8). Specifically, ZOTU_2654 (genus Sporanaerobacter) and ZOTU_52 (genus Tepidanaerobacter) were linked to MAG1 and MAG2, respectively. ZOTU_4243 and ZOTU_6705 were associated with MAG3, classified as Limosilactobacillus reuteri. ZOTU_756, ZOTU_2098 and ZOTU_5627 were mapped to MAG4 and classified as Ruminococcaceae JAAZLM01. All matched MAGs were of high quality, with completeness ranging from 80.24 to 95.28% and contamination ranging from 0.67 to 3.52% (Table S9). These key species were likely to play essential roles in the metabolism of substances during the dynamic process, as they each generally possessed a relatively complete fatty acid biosynthesis pathway and exhibited diverse capacities in carbohydrate metabolism and fermentation (Figs. 4D and S10). For example, MAG1 showed proficiency in butyrate metabolism, while MAG2 contained the complete pathway from glucose to pyruvate to acetyl-CoA to acetate (Fig. 4D). Moreover, both MAG3 and MAG4 matched with multiple key species and demonstrated the ability to metabolize macromolecular carbohydrates such as starch and the ability to produce alcohol. Especially, MAG4 was a generalist in carbohydrate metabolism with the capability to utilize 15 different substrates, including various polysaccharides such as xylan, xyloglucan, cellulose, and others. Therefore, several vital microorganisms organized the system’s metabolic framework through performing crucial functions within the system and possessing a wide range of metabolic capabilities. In addition, the key species directly linked to CH₄ in the network did not show methanogenic functions. Since the network was constructed by selecting the strongest correlations and using the RMT method, the threshold was set at 0.819. Although the correlations were significant (p < 0.05), the links between methanogens and CH₄ were not retained in the network, as their strongest correlation was 0.656.

Discussion

Understanding the strategies of microbial anaerobic meta-metabolism can guide the utilization and regulation of material and energy cycles in various natural and artificial ecosystems, as well as contribute to the maintenance of ecosystem health [1, 2, 4, 8]. However, the complexity of metabolic activities and flexibility in metabolic strategies under different environmental conditions present challenges in elucidating the link between microbial metabolisms and metabolite turnover [20]. Here, we maintained an anaerobic digestion system as a model system to elucidate microbial metabolic rules by monitoring microbial and metabolite dynamics. Instead of being limited to system-wide metabolic performance [64, 65], we adopted a meta-metabolism perspective and focused on both overall performance and individual strategies within the AD system over time. By drawing parallels between metabolite composition and microbial communities, we have investigated the diversity, turnover, community ecology, and bipartite networks of metabolites. However, it is crucial to note that metabolite composition and microbial communities are not entirely equivalent. For example, metabolites do not reproduce or migrate actively, and their interconversion is distinct from microbial processes like reproduction, competition, or extinction [18]. Our method leveraged the approximate characteristics to enhance our understanding of microbial metabolism. Additionally, we highlight the role of microbial function and metabolic traits in the performance of metabolic activities, hence uncovering the underlying mechanisms of co-turnover between microbial community and DOM metabolic composition. Consequently, we proposed three general dynamic rules (Fig. 5) of microbial anaerobic meta-metabolism guiding the DOM conversion. By discovering the intricate relationships and interactions between microorganisms and DOM metabolites, we gained deeper understandings of the metabolic strategies employed in specific systems, facilitating further prediction and regulation efforts.

Fig. 5

The dynamic rules of microbial anaerobic metabolism within anaerobic digestion systems. Three general dynamic rules were concluded as follows: (I) Pathway of deterministic process from operational conditions to microorganisms to metabolites, (II) microorganisms respond to environmental conditions through metabolic trade-offs at community level, (III) microorganisms cooperate based on metabolic complementarities and substance exchange

Pathway of deterministic process from operational conditions to microorganisms to metabolites

The compositions of DOM metabolites in the AD systems were decomposed or synthesized rapidly under a deterministic pattern. The complexity of molecules is generally associated with the overall persistence of organic composition [66, 67]. A recent study on the decomposition of soil plant residue observed a divergent-convergent pattern of microbial meta-metabolism, where complexity increased from 0.5 to 3 years and decreased from 3 to 9 years [16]. In the AD system, we observed a similar trajectory of metabolite diversity, with a divergence-convergence pattern occurring as early as day 3 (Fig. 1C). The rapid biological metabolism in the AD system [25, 26] significantly reduced the time required to observe metabolic patterns in natural environments [16, 17]. Furthermore, the composition of DOM metabolites in parallel bioreactors underwent continuous turnover in a similar trajectory (Figs. 1C, D, E, F and S2). This pattern of constant turnover, previously observed in microbial community succession, suggested that external conditions determine compositional variation [68, 69]. Through assembly process analysis, we demonstrated that metabolite composition was primarily structured by deterministic processes (Fig. 3A). This implied that carbon turnover during the AD process is regulatable and predictable, providing opportunities for engineering interventions in carbon activities through adjustments in process parameters or the addition of microbial agents to modify the microbial community structure or functional expression. However, the deterministic process of metabolite turnover arises from a combination of biotic and abiotic factors, and the relative contributions of these factors are subject to debate [1, 70]. Metabolites are generated from microbial activities [3] but are influenced by complex interactions with local environmental conditions [9, 70]. Consequently, the associations between metabolites’ and microbial communities’ compositions can vary across habitats [22]. Our study provided some evidence for the co-turnover of both microorganisms and metabolites within an anaerobic digestion system (Table S5 and Fig. 3B). Importantly, our results suggested that this co-turnover with microorganisms might be more significant for metabolites than the influences of environmental conditions (Fig. 3B). In addition, a noteworthy observation was the gradual accumulation of stochasticity within the system and a subsequent decline in selection (Fig. 3A). As the selection pressure from the environment increased over time (Fig. S3), this shift from determinism to stochasticity was primarily driven by biotic factors. This shift signified a reduction in microbial metabolic capacity and activity. Consequently, the rapid conversion of molecules (Fig. 1C, D, E, F) and the deterministic processes (Fig. 3A) observed in AD systems were predominantly the outcomes of microbial activity. Meanwhile, microorganisms exhibited sensitivity to environmental pressures and variations (Fig. 3B), implying that environmental factors may regulate the composition and turnover of metabolites by influencing microbial communities. Therefore, microorganisms and DOM metabolites seemed to synchronize their responses to environmental conditions, resulting in coordinated assembly processes (Table S7). As microorganisms responded to dissimilar environmental conditions, the metabolic activities were altered responsively, which in turn affected the composition of metabolites, indicating the environmental deterministic effect within the system was transmitted from the microbial community to the metabolite composition.

Microorganisms respond to environmental conditions through metabolic trade-offs at community level

The overall performance was the consequence of individual metabolic behaviors, and the metabolic patterns exhibited variability across different categories and conditions. In the AD system, the overall performance was characterized by a biodegradation process leading to a reduction in molecular weights (Table 1). A previous study suggested that catabolism drives diversity, while synthesis promotes monotony [3]. However, in our system, continued decomposition and reduced diversity were highly coordinated (Fig. 1C and Table 1). This result indicated that it was impossible to summarize microbial community metabolic patterns in terms of synthesis and catabolism dichotomies, and the metabolic performance was depended on the prevailing conditions and concomitant metabolic trade-off [71, 72]. Meantime, accurately capturing individual metabolic patterns within the context of complex macrometabolic systems remains challenging, at least for the time being [73]. However, the coarse-grained metabolic patterns at the community level were traceable. For instance, lipid-like metabolites followed a distinct trajectory from other categories, with an initial period of biosynthesis (Fig. 1E and Table S4). Microorganisms tend to accumulate lipids as cellular raw materials and energy reserves [74]. However, the efficiency of material transport and energy conversion limits the microbial choice for carbon chain elongation [75, 76], as reflected in the accumulation of C18 substances (Fig. 1F). Therefore, these C18 substances have the potential to serve as an indicator of biosynthesis during the AD process. Furthermore, the degradation of C18 substances occurred in the later stage of the experiment (Fig. 1F), suggesting a shift in metabolic strategy to sustain community functions as resources became depleted (Fig. S3) and thermodynamic conditions changed (Fig. 2D). In addition, the biotransformations performed by microorganisms represented the metabolic strategies [49, 52]. Hence, the temporal specificity of biotransformation (Fig. 2B) can be a specific manifestation of the flexibility in shifting metabolic strategies. The preferential transformation of N-containing substances and the delayed degradation of S-containing substances (Fig. 2B and C) indicated the strategic prioritization of resource utilization by microbial community. Thus, the behaviors of microorganisms and metabolites varied over categories and the environmental conditions by favorable metabolic trade-offs in communities, where the community biotransformation thermodynamics serve as important indicators.

Microorganisms cooperate based on metabolic complementarities and substance exchange

The metabolic capabilities of a microorganism are defined by its genome, encompassing the range of direct metabolic activities it can perform [48, 77]. During the AD process, key microorganisms were identified with capabilities to elongate fatty acids, metabolize various carbohydrate substrates, and produce short-chain fatty acids and ethanol (Fig. 4D). Interestingly, after retaining only the strongest correlations, these remaining key microbes might not be directly involved in the production of CH₄ [78], a pivotal metabolic activity in the system (Fig. 4A). Instead, the modulation of methanogenic metabolism could be influenced by these key microbes, which are more proactive in responding to environmental conditions [26, 79, 80]. This suggested that the abundance of methanogens might not be the sole determinant. As previously reported, microorganisms involved in hydrolysis and fermentation exhibited higher functional redundancy than methanogens [26]. Their metabolic activities during the initial steps of AD may constrain subsequent processes through mutual functional dependencies and interactions [26]. This observation opens up the possibility of engineering the functional microorganism structure to optimize methanogenic conditions, showcasing the potential for targeted regulation of carbon flow in AD systems [81]. Metabolic exchanges are ubiquitous in microbial communities [61], and we connected metabolites traits with metabolic exchanges. During the experiment, microorganisms selected organic matters with specific molecular characteristics as manifested in aggregating molecules with distinct traits by diverse functional microorganisms within the bipartite network (Fig. 4A and B). These molecular traits encompass elemental compositions, bioavailabilities, and barriers to conversion, all directly correlated to molecular behaviors in the process of synthesis, degradation, and diffusion. Consequently, the modular arrangement of molecules with varying characteristics around functional microorganisms (Fig. 4B) illustrated the dependence of microbial metabolic functions on molecular traits (Fig. 3C). Furthermore, all network connectors between microbially dominated modules were contributed by molecules (Figs. 4A and S9). Functional microorganisms likely rely on these molecules for substance exchange, hence achieving metabolic complementarity among these microorganisms and functional complementarity among the network modules. Meantime, microorganisms preferred the biodegradation of larger molecules in a synergistic pattern (Fig. 4C), leading to their swift depletion (Table 1 and Fig. S5). This preference for degrading larger molecules, which provides more resources and energy, supported these microorganisms’ essential activities while concurrently enhancing the resource availability of DOM metabolites (Table 1) for mutual benefit with other microorganisms [29, 82, 83]. Leveraging the shared material pool within the system [78], microorganisms cooperated and influenced each other through metabolic selection, culminating in the formation of an active interaction network. Due to substance exchange and modifications of physical conditions, microorganisms realized synergistic cooperation during anaerobic metabolic activities, driving the rapid biodegradation of larger molecules, the government of methanogenesis, and the maintenance of system harmony.