Introduction

Microbial communities drive all biogeochemical processes on Earth through spatiotemporal resource partitioning among co-existing members or ecotypes that are ecologically differentiated by their preference for a particular environmental context like seasonal fluctuations or physical attributes like free-floating (planktonic) or attached lifestyles1,2. Ecotype differentiation, which leads to resource partitioning, has long been known to be a key mechanism by which microbes divide and conquer diverse niches. There are numerous examples of how even within a species (e.g., species of Prochlorococcus3, Vibrio2, etc.) resources are partitioned in a single habitat through ecotype differentiation. Each ecotype is specialized to maximally utilize some subset of resources, each occupying a different niche. In fact, only recently using comparative genomics have we gained a better understanding of how finely niche space is divided, even in what might appear to be a homogeneous environment. This adaptive diversification is the underpinnings of the remarkable diversity of microbial life now being more revealed through the developing science of metagenomics. It has recently been estimated that 20–80% of the cells in subsurface environments exist as biofilms exhibiting, relative to their planktonic counterparts, distinct physiologies and phenotypes, such as increased resistance to antimicrobials, heavy metals, desiccation and substrate deprivation4,5,6,7. The local genotype interactions among these ecologically differentiated strains confer emergent properties such as stability and cooperativity to the microbial community8. While there is evidence that extended periods of evolution may eventually drive sympatric speciation of ecologically differentiated strains9,10,11, there is scant evidence of early events, mechanisms, and time frames over which ecological differentiation manifests in a microbial community of clonal isolates of two or more organisms.

We have discovered early events in ecological differentiation of a nascent community of sulfate reducing bacteria (SRB) and methanogens engaged in a type of mutualism called syntrophy, which is responsible for transformation of > 1 gigaton/year of C into CH412. This particular type of syntrophy can be either obligate or facultative, depending on the availability of sulfate. Sulfate respiration (SR) is energetically favorable compared to syntrophic growth, which requires interaction with a methanogen in a sulfate-depleted and an energy limited anaerobic ecological niche13,14,15,16,17,18. The two-organism synthetic community (SynCom) of a SRB, Desulfovibrio vulgaris Hildenborough (Dv)19, and an archaeal methanogen, Methanococcus maripaludis S2 (Mm)20, has been used as a model to study syntrophic interactions in the laboratory16,17,18,21,22,23,24,25. When cocultured in the absence of sulfate, Dv mainly oxidizes lactate to acetate, minimal amounts of formate and H2, but the build-up of H2 partial pressure eventually becomes growth inhibitory. Consumption of H2 by Mm to produce methane (i.e., hydrogenotrophic methanogenesis) supports its growth, while reducing H2 partial pressure and alleviating Dv growth inhibition, making the overall reaction energetically favorable13,26. Previous studies have indicated that H2 from Dv is the only electron carrier when paired syntrophically with Mm17,18,27. However, flux balance analysis (FBA), using a constraints-based metabolic model of syntrophy, predicted that in the presence of H2 there could be exchange of formate from Dv to Mm, but that by itself formate cannot serve as the sole electron carrier17. Alternatively, alanine could serve as a C and N exchange intermediate from Dv during syntrophic growth with Mm, but this only occurs in the absence of ammonia16,28. Ultimately, the energetically favorable reduction of H2 to methane, by the SRB and methanogen allows these microbes to access a niche previously uninhabitable to both independently, and together drive one of the major processes of methane production on Earth.

Early on in the laboratory evolution of syntrophic interactions between Dv and Mm there were periods of instability and extinction events, but stability emerged within 300 generations significantly improving growth rate and yield24. Subsequent analysis discovered that while each evolved partner had individually contributed to improvements in growth characteristics, pairings of evolved clonal isolates of Dv and Mm from the 1000th generation (1 K) synergistically improved overall growth characteristics. However, maximal improvements in growth rate and yield were observed in minimal assemblages of the syntrophic community that were obtained through end-point dilution of the evolved community from the 1,000th generation2,4. Notably, dilution to extinction had generated distinct sub-populations that had segregated mutated genes with high G-scores (i.e., mutations observed in the same gene across multiple independent evolutionary lines) that co-existed for up to 700 generations of laboratory evolution. In fact, some minimal assemblages had enriched variants that were below the limit of detection in bulk sequencing of the parental 1 K generation line. In the absence of distinct functional interactions, we would have expected to see all mutations from the parental line represented proportionally within each minimal assemblage. Thus, these observations suggested that overall growth improvements observed during laboratory evolution had likely emerged from structured interplay of multiple variants across the two organisms, including low frequency variants that had retained SR capability22,23.

Using custom fluidized bed reactors (FBRs) to simulate an ecologically relevant subsurface environmental context for syntrophic interactions in soil (i.e., attached to sediment and free-floating in groundwater) we investigate how the interplay between variants of vastly different abundance and physiological capabilities contribute to improved growth characteristics of the minimal assemblage of the laboratory evolved community. Specifically, we ask whether diversification of genotypes in one environmental context (i.e., planktonic growth) could give rise to genotypes that might thrive in a separate environmental context29,30 (i.e., attached to sediment) to improve overall community characteristics? We perform longitudinal assessments of changes in genetic diversity and transcriptome responses of Dv and Mm across attached and planktonic phases in the context of overall growth and productivity of the SynCom (i.e., lactate utilization and methane production). Further, by developing contextualized genome scale metabolic network models for biogeographically differentiated syntrophic communities, we investigate the dynamic changes in metabolic flux states across the planktonic and sediment phases. In so doing, we follow the segregation of variants that emerged at the origin of a nascent microbial mutualism into biogeographically interdependent sub-communities, quantifying both the emergence and maintenance of new ecotypes. We observe how mutations in a few genes, across the two organisms from a previously evolved community, mechanistically lead to ecotype differentiation in a simulated subsurface environment, likely driving resource partitioning and improving overall productivity of the mutualistic microbial community.

Results

Adaptation of a laboratory evolved syntrophic community to a simulated sub-surface environment

Custom anaerobic FBRs were designed as an analog for a subsurface environment, wherein sediment and groundwater are in flux and free to exchange metabolites and biota31,32. Three independent FBRs were inoculated concomitantly with a minimal assemblage of the Dv-Mm SynCom obtained through end-point dilution of the mutualistic community adapted to obligate syntrophy conditions through 1000 generations of laboratory evolution22,23,24. Each FBR contained 75 g of sediment (210–297 μm crushed quartz) and approximately 350 mL (FBR and reservoir) of lactate coculture medium (Supplementary Data 1) without sulfate to impose obligate requirement of growth by syntrophy (Table 1 and Fig. 1A). The sediment bed in each FBR was briefly fluidized by recirculating the lactate coculture medium from the bottom of the reactor and up vertically through the column to mix cells throughout the reactor. The SynCom was cultured in batch mode for approximately 48 h, as done previously15. After 2 days of growth, the FBRs were set to continuous fluidization (Day 1–6) (350 mL min−1) with a duty cycle of 1 h on and 1 h off for the remainder of the experiment. The triplicate reactors were operated anaerobically (Dissolved Oxygen: DO ~ 1 mg/L) at pH ~7.2 for approximately 175 h under semi-continuous batch culture conditions with supplementation of fresh lactate medium at 7.8 mL h−1, with a net dilution rate of 0.022 h−1 (Table 1, Supplementary Fig. 1), which was similar to conditions used for growing biofilms (0.017 h−1)15.

Table 1 Operational parameters for anaerobic FBR arrays (n = 3)
Fig. 1: Partitioning of attached and planktonic syntrophic communities in FBRs.
figure 1

A A schematic of the custom FBRs developed to simulate biphasic growth of microbial communities by recirculating growth medium upward through a column of sediment. Temporal growth dynamics of planktonic (B) and sediment attached (C) phases of the microbial community. Metabolite profiles provide confirmation syntrophic growth via the oxidation of lactate to acetate by Dv (D) and the production of methane by Mm (E). The grey opaque background bars indicate the first 48 h of growth in batch culture mode (beginning of Day 1 through the end of Day 2) in FBRs prior to fluidization (i.e., Day 1 through Day 6, indicated in bold typeface); grey triangles represent timepoints (from Day 1 through Day 6), when SynComs from planktonic and attached phases were harvested for DNA and RNA profiling. All data presented as mean values +/− standard deviation (95% Confidence Interval) across 3 replicate FBRs (n = 3). Source data are provided as a Source Data file.

Biomass in the planktonic phase increased at a steady rate with a maximum growth rate of 0.23 d−1 and an average doubling time of ~3 days (Fig. 1B). Using total protein from washed sediment as a proxy for biomass we determined that the attached community accumulated at a slower rate of 0.087 d−1 (doubling time: ~8 days) (Fig. 1C). Syntrophic coupling between Dv and Mm was evident in sustained lactate oxidation to acetate, with the subsequent production of methane (Fig. 1D, E). While significant increase in methane production was observed within 72 h, the levels were steady at ~0.75 mM for the remainder of the experiment. However, it should be noted that since the reactor headspace was allowed to vent through a one-way check valve to prevent pressure from exceeding ~1 atm, the methane measured does not reflect the total amount produced throughout the experiment. Together, these results demonstrated that the coculture established syntrophic growth in the simulated subsurface environment within the FBRs.

Heterogeneity of population structure and genetic variability across growth phases

Total DNA and RNA was sampled daily for 6 days once the reactor fluidization was turned on (Fig. 1. Grey triangles) at approximately 48 h (Fig. 1. grey bar) to evaluate population dynamics and functional heterogeneity of Dv and Mm across attached and planktonic growth phases. Based on normalized numbers of metagenomic sequencing reads corresponding to each organism, we determined that Dv was the dominant partner in both phases but community composition was significantly different between the planktonic and attached communities (Fig. 2A). For instance, Mm in the planktonic phase nearly doubled in relative abundance from day 1 to day 2 but thereafter accounted for about a third of the SynCom composition (i.e., ~3:1 Dv:Mm). In contrast, the relative proportion of Dv in the attached SynCom was ~18 times that of Mm. The relatively lower abundance of Mm in the attached SynCom is likely because Dv attachment occurs first and acts as a scaffold for Mm to join the biofilm, which takes ~25 days for full maturation15. However, it is also likely that variants of both organisms within the evolved SynCom have different propensities of adaptation to syntrophic interactions in either the planktonic or the attached phases. Regardless, these findings suggested different mechanisms of syntrophic interactions had likely manifested across the attached and planktonic phase SynComs.

Fig. 2: Population structure and dynamics between attached and planktonic communities.
figure 2

A Longitudinal change in relative abundance of Dv and Mm, across both phases, calculated from coverage of reads binned by GC content across all samples (n = 35). Error bars represent mean + standard deviation (95% Confidence Interval). B Frequency distribution of Dv mutations (> 0%, < 100%) between planktonic (mean 72.3%, dashed line) and attached populations (mean 79.9%, dashed line); asterisk denotes significant difference (p-value = 1.3e-17, two-sided t-test). C Distribution of Dv mutation frequencies by day; asterisk denotes significant decrease relative to day 1 (p-value < 0.05; two-sided t-test), specifically on days 3 through 6, while attached communities show a shift to higher frequency mutations. D NMDS plotting for all SynCom samples across phases (n = 35). Arrows show ordination of specific mutations in NMDS space representing the optima habitat (showing top 3 for Dv and Mm). Large circles represent the centroid of the two phases. E Heatmap of mutations clustered (Pearson correlation) across phases, days and replicate reactors. Selected variants are highlighted. Impact Prediction: “high impact” mutations include gain or loss of start and stop codons or frameshift mutations; “moderate impact” mutations include deletion or insertion of codons and nonsynonymous changes in coding sequence; intergenic mutations are classified as “modifier mutations”37. F Significantly decreased expression (two-sided t-test) of MMP1667 (archaellin flaB) over days 3 through 6 in the planktonic phase (n = 11) due to the loss-of-function mutation in archaellum regulator (Mmp1718p.R23fs) compared to the attached phase (Mmp1718wt i.e., no loss-of-function) (n = 11). Boxplots: median center line, box limits are upper and lower quartiles; whiskers, 1.5× interquartile range. G Change in average mutation frequency of nonsynonymous variant in galUP32S (DVU1283); dashed line indicates linear regression and R2 is the correlation statistic for longitudinal trends for frequency change of each phase. The ribbon represents standard error (95% Confidence Interval). H Scanning electron microscopy images show EPS overproduction by evolved Dv with the nonsynonymous galUP32S mutation, but not by ancestral Dv-galUwt cells, which was confirmed across three independent evolved lines. Source data are provided as a Source Data file.

We analyzed the metagenomic data using sequence alignment and variant calling to quantify mutation frequencies of Dv and Mm variants in SynComs across the two phases. The overall distribution of mutation frequencies of Dv was significantly lower in the planktonic as compared to the attached SynCom (Fig. 2B). Further analysis revealed that the mean frequency distributions of Dv mutations were statistically similar across the planktonic (78.5) and attached (76.0) phases on day 1, but over time the mean frequency distribution in the planktonic phase trends progressively lower, while it remained relatively stable in the attached phase (Fig. 2C). This progressive decrease in Dv mutation frequencies over time in the planktonic phase was likely due to wash out of genotypes with growth rates lower than the dilution rate (Supplementary Fig. 2A, B). In contrast, the overall mutation frequencies of the Mm population remained stable in both phases. We performed nonmetric multidimensional scaling (NMDS) analysis to investigate if there was evidence of ecological differentiation33 by specific Dv and Mm variants (Fig. 2D). While most Dv variants, including DVU0168 and DVU1295, oriented towards attached habitats, most Mm variants, including MMP1227 and MMP1718, showed clear ordination towards planktonic growth (Fig. 2D shows the top 3 variants of each Dv and Mm). Strikingly, the MMP1718 variant and was only observed 2 of 18 samples in the sediment phase (Fig. 2E). Consistent with the demonstrated role of MMP1718 as a transcriptional activator of the fla operon34, a loss-of-function frame-shift mutation at arginine 23 in this gene (Mmp1718p.R23fs) resulted in decreased expression of archaellum genes (specifically, flaB2) in the planktonic phase (Fig. 2F), which may have disrupted surface attachment of Mm35,36. While this was also consistent with the observation that the Mm-Mmp1718p.R23fs variant had likely enhanced fitness during laboratory evolution, which was performed through sequential transfers of the planktonic phase evolved cocultures22, it was intriguing that the WT allele remained even after end point dilution and ultimately swept through the Mm population in the attached phase.

Of the total 22 (15 Dv and 7 Mm) mutations, nearly half (n = 10) were predicted by SnpEff to have high impact, meaning they disrupted function of the encoded protein, by introducing a gain or loss of start and stop codons or frameshift mutations37. The other mutations were predicted to have moderate impact, such as through introduction of a nonsynonymous change, meaning deletion or insertion of triplet codon(s)37 (Fig. 2E). One mutation in an intergenic region (IG-2053158) in Dv, was classified as a modifier, given the challenge of predicting its consequence37. Using linear regression across the three independent reactors (Supplementary Fig. 2B), we discovered that while frequencies of most variants were stable in the attached community, mutations in a few Dv genes showed preference for the attached phase. In particular, a nonsynonymous variant that changed a non-polar proline to a polar serine in a UTP-glucose-1-phosphate uridylyltransferase (i.e., galUP32S), increased in frequency in the attached phase (R2 = 0.55) (Fig. 2G). Further investigation, using scanning electron microscopy, demonstrated that the nonsynonymous mutations in galU had likely resulted in overproduction of extracellular polysaccharides (EPS), which may have facilitated attachment to sediment (Fig. 2H)38,39,40. Thus, the NMDS analysis uncovered evidence that specific variants of both organisms had likely led to niche differentiation of the SynCom ecotypes across the planktonic and attached phases.

Distinct mechanisms of adaptation and mutualism manifest across co-existing planktonic and attached syntrophic communities

We analyzed longitudinal transcriptome changes to understand adaptive consequences of ecological differentiation on the nature of syntrophic interactions across the planktonic and attached communities. A similar proportion of genes in Dv (536 genes, 14% of total genes), as compared to Mm (295 genes, 17% of total genes), were differentially expressed over time and across the two phases. There was a cascading pattern of changes in gene expression across the planktonic phase cells of Dv, but the transcription of these genes was relatively stable with sustained up- or downregulation in the attached phase cells (Fig. 3A). These expression patterns of Dv genes likely reflect adaptation of a stable, slow growing, sessile community to the attached phase, as compared to faster growing planktonic cells that experience greater environmental dynamics. In stark contrast, although Mm transcriptional response was relatively uniform over time in both phases, there was a significant shift in expression levels of 173 genes during transition from day 1 to day 2 (Fig. 3A), which is consistent with how the two organisms establish syntrophy, especially when attaching to surfaces, with SRB colonizing first and laying the foundation for subsequent attachment of the methanogen15.

Fig. 3: Transcriptome profiling shows distinct expression responses of co-existing planktonic and attached syntrophic communities.
figure 3

A Heatmap (rows: genes) of z-scored expression (TPMs) of sequentially ordered significantly expressed genes (Fold Change > 2, p-value < 0.01, DESeq2 Wald test) for each day (planktonic vs attached) and split between up- or downregulated genes for Mm and Dv. All p-values can be found in supplementary data file 3. Right side annotation represents the log10 mean expression level for each gene across all samples (n = 35). B PCA of transcriptomic data of the evolved community from this study and cocultured planktonic WT Dv and Mm cells from a previously published study showing transcriptional states changes resulting from 1000 generations of syntrophic growth and between phases. C Genes encoding central functions (lactate permeases: n = 3, Flagella: n = 13, Dv Hydrogenases: n = 12, Mm High-Affinity Hydrogenases: n = 23, Mm Low-Affinity Hydrogenases: n = 4) for syntrophic interactions between Dv (orange) and Mm (gray) have distinct expression patterns across planktonic and the sediment attached phase. All error bars indicate the standard deviation (95% Confidence Interval) from the mean (marker) across 3 replicate FBRs. Source data are provided as a Source Data file.

We investigated whether the transcriptome state of the SynCom had changed significantly over 1 K generations of laboratory evolution, and after adaptation of the evolved community to attached and planktonic phases in the FBR. Specifically, we performed principal component analysis (PCA) to compare transcriptome states of Dv and Mm before evolution (i.e., planktonic coculture of WT isolates21), after 1 K generations of laboratory evolution (day 1 of culturing a minimal assemblage from the 1 K generation in the FBR), and across the attached and planktonic phases over the subsequent five days of culturing. The PCA showed that transcriptome state of the wild-type SynCom (i.e., before evolution) had changed significantly after 1 K generations of evolution. There was further divergence of transcriptome states of the evolved SynCom upon adaptation to planktonic and attached phases (Fig. 3B).

Expression patterns of key genes for motility and attachment, lactate utilization, and methanogenesis revealed distinct physiological states of the two organisms across the two phases. For instance, expression of 13 flagellar biosynthesis genes in Dv were constitutively down- or upregulated in the planktonic and attached phases, respectively, which was consistent with the known role of these genes in driving attachment35,36 and partner selection41,42. A similar expression pattern was observed for Mm, although in this case, the constitutive downregulation of the archaellum encoding genes could be attributed to selection of loss-of-function mutations in the cognate transcriptional activator (i.e., Mmp1718p.R23fs) in the planktonic phase. Further, all three lactate permeases in Dv were upregulated in the planktonic phase, especially day 4 onwards when there was an accumulation of biomass and a decline of available lactate for oxidation, which is the primary mechanism for generating H2 as a reductant for establishing syntrophic growth. Interestingly, upregulation of DVU2451, which can functionally compensate43 for the primary lactate permeases DVU2110 and DVU302644, was likely a lactate scavenging mechanism by the planktonic Dv cells. By contrast, the attached Dv cells appeared to be not lactate limited as all permeases were constitutively downregulated in sediment.

Some methanogenesis genes in Mm were upregulated in the planktonic phase (e.g., formylmethanofuran dehydrogenase and F420-reducing hydrogenase), whereas a different set was upregulated in the attached phase (e.g., H2 methylenetetrahydromethanopterin dehydrogenase). Notably, while the low affinity hydrogenases45 in Mm were upregulated in the attached phase, the high affinity hydrogenases were upregulated in the planktonic phase. This finding is consistent with the expectation of significantly higher H2 partial pressure in the microenvironment surrounding Mm in the attached phase, due to physical proximity to Dv cells, which were also in relatively higher abundance (Dv:Mm ~18:1) in sediment. Conversely, lower cell density and relative abundance of Dv likely generated lower H2 partial pressure in the planktonic phase, triggering upregulation of the high affinity isoform of the hydrogenase in Mm as a mechanism for H2 scavenging (Fig. 3C). Thus, results of the expression analysis illustrated how ecological differentiation driven by differential segregation of mutants (e.g., Dv-galUP32S and Mm-Mmp1718p.R23fs) led to further niche partitioning through differential expression of key genes across the two phases of growth within the same reactor, while still maintaining nutritional cooperation (i.e., syntrophy).

Flux balance analysis using a curated genome-scale metabolic model for syntrophy reveals differential reaction fluxes across the two phases

A genome-scale metabolic network model of syntrophy was extended and curated to investigate reaction flux through the two members of the syntrophic community. In brief, metabolic models iJF477 (Dv)46 and iMR539 (Mm)47 were updated by curating 1014 and 688 reactions, respectively, and subsequently integrated into a genome-scale model of syntrophic interactions by setting interfacial constraints for interspecies exchange of 54 metabolites, including formate, H2, and alanine16,18. Further, the new syntrophy model (iSI1283) was contextualized to planktonic and attached phases and across phases by using the GIMME algorithm48 to exclude reactions encoded by genes that were expressed below a threshold (Fig. 4A, Table 2).

Fig. 4: An integrated multi-species metabolic model reveals distinct metabolic flux states across attached and planktonic phases during syntrophic growth.
figure 4

A Schematic of a constraints-based metabolic model for syntrophy (iSI1283) obtained through the integration of genome-scale metabolic networks of Dv (iJF477) and Mm (iMR539). The syntrophic model highlights the general exchange of key metabolites for hydrogenotrophic methanogenesis and the required in vivo exchange, and transport reactions between intracellular and extracellular compartments. B In silico prediction of relative biomass of Dv and Mm before (small, transparent points, with dashed purple line) and after (large, solid points, with solid black line) constraining Dv and Mm biomass ratios. The ribbon on each linear regression (lines) represents standard error (95% Confidence Interval). C t-SNE analysis of flux states of Dv and Mm over time and across planktonic and attached phases. Relative fluxes towards biomass of each organism, predicted by FBA using the syntrophic model with methane production as the objective function (D, E). Model predicted rates of key metabolites associated with methanogenesis for each phase and day: (F) lactate uptake by Dv, (G) H2 production by Dv (Total bars: hashed and non-hashed), H2 utilization by Mm (non-hashed), and (H) amount of methane produced by Mm. I A heuristic model for the predicted differences in relative flux (mMolgDCW−1h−1) between planktonic and attached cells and ecotypes (Dv-galUwt —Dark orange: Mmp1718p.R23fs —Light grey and no archaellum; Dv-galUP32S —Light orange and thick membrane: Mmp1718wt —grey with archaellum). Arrow weights reflect the amount of relative flux predicted through or to each organism and phase. Source data are provided as a Source Data file.

Table 2 Summary of reactions, genes, and metabolites in metabolic models for Dv and Mm, individually, and together in syntrophic interaction

Relative abundance of the two organisms in the attached and planktonic communities was accurately predicted (R2 = 0.77) by the unconstrained model when it was contextualized with gene expression data, with upper bounds on total biomass of Mm, demonstrating that it had captured general trends of metabolic interplay and community structure within and across phases (Fig. 4B). For further investigation of these evolved strategies, the experimentally quantified relative abundances of both Dv and Mm (Fig. 2A) were included as constraints for prediction of syntrophic growth of both organisms, which significantly improved model accuracy (R2 = 0.92 from 0.77) (Fig. 4B, Supplementary Fig. 3). Further, t-SNE (Stochastic Neighbor Embedding) analysis using biomass constrained syntrophy models demonstrated that while flux states of the SynCom were indistinguishable on day 1, they became increasingly divergent over time as the communities partitioned across the two phases (Fig. 4C). We used the iSI1283 models to uncover relative contributions of SynComs in the attached and planktonic phases (i.e., the relative flux values given equal amounts of dry cell weight— DCW) to overall productivity of hydrogenotrophic methanogenesis in the FBRs. Notably, when optimizing for methane production (i.e., syntrophic growth) the relative flux towards Dv biomass was higher in the attached phase for all days (Fig. 4D). In contrast, the relative flux towards Mm biomass was higher in the planktonic phase (Fig. 4E). This analysis predicted that the relative rate of lactate utilization by Dv was on average significantly higher in the attached phase compared to the planktonic cells (p-value = 0.034; t-test) (Fig. 4F). The production rate of H2 (Dv) was stoichiometrically proportional to rate of lactate utilization. Consistent with findings reported in previous studies, the model also predicted that H2 was the dominant electron exchange intermediate for both planktonic and attached Mm, whereas electron transport through exchange of formate was near zero (Supplementary Fig. 4)17,18,27. Interestingly, the contextualized model predicted that a significantly higher fraction of H2 produced in the attached phase was unutilized (~59.4%), as compared to the planktonic phase (15.1%) (p-value = 5.3e-5; t-test; Fig. 4G). Methane production was significantly higher in planktonic phase, relative to the attached phase, likely due to the excess H2 from both phases and the higher relative proportion of Mm (Fig. 4H). The FBA prediction of low efficiency of H2 utilization in the attached community was consistent with the higher expression level of low affinity hydrogenases in Mm. The excess H2 release by the attached community was likely scavenged by high affinity hydrogenases towards significantly higher rate of methanogenesis by planktonic Mm (Fig. 4I).

Discussion

Our study demonstrates that, even within a simple two organism community that originated from clonal isolates, ecotypes adapted to different biogeographies can arise very rapidly (300–1000 generations), introducing complex structure and division of labor, thereby increasing the overall productivity of the community. In the soil subsurface, microbial communities can be more productive in attached (sediment) relative to planktonic phase35,36,49. Two explanations have been put forth to explain this phenomenon: first, that the particulate material in sediment is largely organic and nutrient rich, or inorganic material with adsorbed nutrients36,50,51. Second, it has been argued that an attached community in a flowing environment has greater access to soluble nutrients due to change in the depth of the diffusive boundary layer - the greater the flow, the thinner the layer and greater the flux of nutrients through the boundary layer35,52,53,54. Both of these mechanisms may be at play in the FBRs and responsible for higher relative rate of lactate utilization by Dv in sediment, also explaining why the Dv-galUP32S mutant might have benefitted from attachment. By contrast, Dv cells in the planktonic phase likely went with the flow and were therefore relatively less efficient at accessing and utilizing lactate, although they leveraged their diverse physiological capabilities as generalists by continually adapting to changing availability of nutritional sources as reflected in the cascade of changes in expression of large numbers of genes. The high amounts of H2 produced by Dv generated a higher local H2 partial pressure in sediment, however H2 was not completely utilized by attached Mm cells despite their physical proximity. Attached Mm cells upregulated low affinity hydrogenases to aid H2 consumption and reduce H2 partial pressure in the microenvironment of attached cells to prevent feedback growth inhibition of Dv. Interestingly, in the attached phase, Mm was also less efficient at transforming H2 to methane, and as a result excess H2 was released into the fluid environment, where planktonic Mm used high affinity hydrogenases to capture/scavenge H2 and use it towards copious amounts of methane production. Indeed, Mm existed in higher relative proportion to Dv in the planktonic phase as compared to sediment, and this also explains why the loss-of-function mutation in the archaellum regulator (Mmp1718 p.R23fs) conferred higher fitness in the planktonic phase (Fig. 4I). The enrichment of certain ecotypes demonstrated their adaptation to the sediment phase through attachment and to the planktonic phase by doubling at a rate (~3 days) that exceeded the dilution rate (~3.7 days). In other words, despite the mixing of genotypes; the ecotypes that had differentially adapted to sediment and planktonic phases through attachment and fast growth, respectively, segregated into distinct communities that were functionally interdependent.

It is clear that interplay of ecotypes of the two organisms across the planktonic and attached phases greatly improved cooperativity and overall productivity (i.e., biomass and methane production) of the entire mutualistic community22. But the emergence of a sediment adapted evolved ecotype in this interplay was perplexing because the laboratory evolution was specifically designed for selection of mutations that optimized syntrophy in a planktonic growth state through serial transfers of liquid coculture. In other words, it is intriguing how mutations that improved characteristics of the attached community arose and how they remained at low frequency. For instance, it is understandable that the Mmp1718p.R23fs mutation was selected because it promoted better planktonic growth of the methanogen, but it’s puzzling how the ancestral allele persisted in the population at low frequency, even after end-point dilution. Our results show that the ancestral Mmp1718wt allele may have benefited from improved physical interactions with Dv, especially upon appearance of the EPS overproducing galUP32S mutant. The Dv-galUP32S mutant may have segregated with Mm-Mmp1718wt as attached communities on the walls of Balch tubes, which begs the question of how they persisted despite multiple serial transfers of the liquid culture, especially given that establishment of syntrophy in a biofilm takes significantly longer than the period of culturing between transfers, and the growth rate of cells is significantly slower in an attached state than in a planktonic state15. Perhaps, the Dv-galUP32S and Mm-Mmp1718wt variants survived as free-floating cell-to-cell interactions (i.e., attached partners)41. Regardless of how these traits persisted through serial transfers of the planktonic phase, our findings suggest that even in well mixed environments, cell aggregation through biofilm formation may improve emergent properties by providing spatial structure to the community55.

Our findings suggest that the planktonic and attached lifestyles of the two organisms in a syntrophic association are intricately intertwined such that optimization of one lifestyle drives co-optimization of the other. This hypothesis is supported by our findings that sediment attachment of the Dv-galUP32S mutant was associated with improved lactate utilization, and may have required physical interactions with Mm-Mmp1718wt for reducing H2 partial pressure in its microenvironment, and conversely, improved planktonic growth of the Mm-Mmp1718p.R23fs mutant was associated with more efficient H2 utilization and also likely benefitted from metabolite exchange with free Dv cells that continually adapted to changing conditions of the flowing medium. Nevertheless, the interplay of physiologies associated with planktonic and attached ecotypes of the two organisms is likely ingrained in the molecular networks of the two organisms, such that mutations affecting one physiological state potentiated selection of mutations in genes associated with the other linked physiological state. Future experiments, could provide important parameters, such as diffusion constants and uptake rates for each extracellular metabolite that takes into account the complex hydrodynamics across gradients and boundary layers, relative abundance of Dv and Mm, and thickness of the attached community. These experimentally determined values will be critical for parametrizing and modeling complex cross-feeding interactions between the planktonic and attached communities of biogeographically differentiated ecotypes56.

Our results demonstrate that physical segregation of specific genotypes that arose over short evolutionary timescales (i.e., 300–1000 generations) led to ecotype differentiation and resource partitioning wherein the SynComs manifested distinct physiological states and vastly different community structures across the attached and planktonic phases. Yet there was meaningful interplay between the attached and planktonic communities that likely contributed to overall higher productivity. Thus, our findings shed insight into the evolutionary mechanisms by which microbial communities leverage ecotype differentiation to divide and conquer spatiotemporally complex environments. Remarkably, our findings also show that ecotype differentiation occurred within 300–1000 generations by mutations in as few as 15 genes in Dv and 7 in Mm, which is consistent with observations in natural microbial communities in the environment11,57,58. Importantly, ecotype differentiating mutations in the SRB and the methanogen, respectively, had previously calculated high G-scores during laboratory evolution experiments, which is a measure of reproducibility of selecting mutations in the same gene across many independent evolution lines22. Evidence of high G-scores indicates that the ecotype differentiating mutations (including the Mmp1718p.R23fs and the galUP32S) were observed independently and at high frequency across multiple evolution lines demonstrating further that this mechanism of ecotype differentiation was fairly common and also that the attached community (primarily Mm-Mmp1718wt and Dv-galUP32S) had co-optimized with the planktonic community across multiple lines. We predict that prolonged laboratory evolution of such a microbial community with co-existing but physically segregated genotypes that exclude each other across habitats could eventually lead to sympatric speciation2,10,11,59.

Methods

FBR setup, operation, and inoculation with SynCom

We designed custom glass anaerobic fluidized bed reactors that operate by recirculating media from the bottom of the reactor to the top creating a vertical velocity through the column to create a bed of fluidized sediment. The anaerobic FBRs were arrayed in triplicate in a temperature-controlled room (~30 °C) and filled with 75 g of 210–297 μm crushed quartz (Sigma-Aldrich: 50–70 mesh Quartz) to act as the sediment bed for microbial attachment, which sits above a diffuser plate (P0 Frit: 160–250 μm) allowing only the flow of liquid generating sediment fluidization. All three reactors were operated at the same time in a temperature controlled room held at ~30 °C. Reactors were filled with sterile lactate coculture medium (Supplementary Data 1) (no sulfate and Resazurin anaerobic indicator) to an operating volume of approximately 350 mL. We then used flow meters (~200 mL min-1) to flush each reactor with an 80%:20% N2:CO2 gas mix through 0.2 micron filters for approximately 15 min to ensure anoxic conditions prior to inoculation. Peristaltic pumps were used for recirculation to generate fluidization at 350 mL min−1 with a duty cycle of 1 h on 1 h off. Fresh lactate coculture medium was added (0.13 mL min−1) to each reactor at the same duty cycle as the rate of fluidization. The reactors were inoculated with a minimal assemblage of Dv and Mm that was obtained from a previous study22,23,24 through end-point dilution of a laboratory evolved coculture after 1000 generations of obligate syntrophic growth. Specifically, a glycerol stock of the minimal assemblage was revived in anaerobic Balch tubes containing lactate coculture medium and subjected to three successive transfers to reestablish growth and then transferred to a 60 mL serum bottle for inoculation into FBRs. A 15 mL aliquot of optical density (OD600) 0.3 culture from a single serum bottle was used to inoculate each FBR containing anaerobic lactate coculture medium. FBRs were briefly fluidized to mix cells throughout the reactor and allowed to grow in batch mode for approximately 48 h to allow biomass accumulation before switching to fluidization mode (Days 1–6) for the rest of the experiment.

Sample collection and analysis

While reactors were fluidizing, approximately 15 mL of planktonic culture was sampled anaerobically using a sterile 20 mL syringe and needle. Planktonic fractions were split for OD600 analysis (1 mL), protein extractions (2 mL), nucleic acid extractions (4 mL), nutrient analysis, pH (2 mL), and dissolved oxygen (DO) (2 mL). Protein and nucleic acid samples were kept on ice until pelleted. In microcentrifuge tubes, 2 mL of planktonic culture was spun down and pelleted at 12,000 × g for 10 min at 4 °C. The decanted supernatant was saved (−20 °C) for organic acid (lactate and acetate) analyses, while the cell pellets for protein and nucleic acids were flash frozen in liquid N2 and stored at −80 °C until extraction. Lactate and acetate samples were diluted with milliQ H2O and measured using ion chromatography60. A summary of monitored parameters over the experimental time-course is provided as Supplementary Table 1.

To sample the communities attached to sediment, we created sample chambers made from Balch tubes that were filled with 100% N2. While reactors were fluidizing, sediment was sampled anaerobically by pulling a vacuum from a sampling port of the FBR into a sterile anaerobic Balch tube that was on ice. Immediately following sediment sampling from the FBR, samples (on ice) were taken directly to a 4 °C cold room for washing. The sediment samples were transferred to 15 mL conical tubes and excess planktonic culture was decanted into a waste container. Approximately, 10 mL of pre-chilled PBS (4 °C) was added to the conical tubes containing sediment was washed for 5 min on a benchtop tube rotator. The PBS wash solution was then decanted from the sediment samples and the 5 min wash step was repeated with fresh chilled PBS. After the second wash, sediment samples were split into three 2 mL microcentrifuge tubes (Protein, DNA, and RNA) and any excess wash solution was pipetted out of the samples were immediately frozen in liquid N2.

We prepared glass vials filled with 100% N2 and crimped closed with rubber septum and aluminum seal to store sampled headspace from the FBRs to be analyzed via gas chromatography (GC). Immediately prior to sampling reactor headspace, 2 mL of 100% N2 was removed from the ethanol washed storage vials using a small diameter needle and 2 mL syringe. The same needle and syringe was then used to sample 2 mL of reactor headspace through the ethanol sterilized headspace sampling port and reinjected into the storage glass vials for quantification of methane using GC. Methane was measured using an SRI 8610 C GC (SRI Instruments, Torrance, CA) equipped with a 6’ silica gel column (SRI Instruments, Torrance, CA) at an oven temperature of 80 °C and a flame ionization detector operated at 385 °C. The carrier gas consisted of >99.998% N2 gas (Praxair, Danbury, CT) at 20 mL min−1, >99.5% H2 gas (Praxair, Danbury, CT) at 25 mL min−1, and air supplied via an internal pump at 250 mL min−1.

DNA, RNA, and protein extractions

Nucleic acids were extracted from planktonic pellets and sediment samples using the MasterPure Complete DNA and RNA Purification Kit (Epicentre: #MC85200). Nucleic acid samples were then split between DNA and RNA samples. DNA samples were treated with RNase and RNA samples were treated with DNase. All samples were quantified using the Qubit DNA or RNA High Sensitivity (HS) kits (Invitrogen). Protein extractions for planktonic cell pellets were performed according to Thermo Scientific B-PER Complete Bacterial Protein Extraction Reagent protocol. For sediment samples, sediment was transferred to a 1.5 mL eppendorf tube and incubated with B-PER Bacterial Protein Extraction Reagent and rotated for 15 min to extract protein. After protein extraction, sediment was transferred to a glass petri dish and dried in a 70 °C oven overnight to quantify the mass of the sediment. All protein samples were quantified using the Qubit Protein Broad Range Assay and normalized to dry sediment weight.

DNA and RNA library preparation and sequencing

DNA sequencing libraries were prepared using the Nextera XT DNA library prep kit, while RNA samples were rRNA depleted using RiboZero Plus kit followed by library construction using the Truseq Stranded Total RNA library prep kit. Both DNA and RNA libraries were sequenced using a mid-output 2 × 75 bp kit for 150 bp fragments using the Illumina NextSeq 500 platform. MultiQC61 output from DNA sequencing was used to quantify relative abundance from metagenomes of the SynCom based on percentage of reads per GC content for Dv (~63.2%) and Mm (~33.1 %) for each sample (n = 35).

Identification of mutations

In this study, we used the same sequence alignment and variant calling pipeline generated in our previous work by ref. 22 (see also https://github.com/sturkarslan/evolution-of-syntrophy62). This pipeline included quality control and trimming of the raw sequencing reads in fastq format by using Trim Galore software63 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore). The alignment of the quality trimmed sequences to reference D. vulgaris (Genbank assembly: GCA_000195755.1.30 [https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_000195755.1/]) and M. maripaludis (Genbank assembly: GCA_000011585.1 [https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000011585.1/]) genomes and subsequent processing steps before calling the variants was done by following The Genome Analysis Toolkit (GATK)64 best practices. Briefly, reads were first aligned to the reference genome using Burrows-Wheeler Alignment Tool (bwa)65 (version 0.7.17-r1188) in paired-end mode. The resulting alignment files in the SAM format were converted to BAM files, sorted and indexed by using Samtools version 1.966. BAM files were marked for duplicates using Picard Tools (http://broadinstitute.github.io/picard/) (version 1.139), and local realignment around indels was performed to identify the most consistent placement of reads relative to the indels. Variant calling was performed independently by using three different algorithms including GATK UnifiedGenotyper, Varscan 267 (version 2.3.9) and bcftools from Samtools package. The default parameters were used for UnifiedGenotyper, whereas for Varscan parameters were --min-coverage 8 --min-reads2 2 --minavg-qual 30 and bcftools parameters were -vmO -s LOWQUAL -i’%QUAL > 30. A variant was included in the analysis only if it is simultaneously called by at least two of the callers. The resulting variants were annotated using SnpEff tools37 (version 4.3). Variant analysis, NMDS, and linear regressions were done using custom scripts and a list of compiled variants is described in Supplementary Data 2.

Calculation of G-scores

In the current study, we refer to G-scores (goodness-of-fit) for ecotype differentiating mutations (including the Mmp1718p.R23fs and the galUP32S). The G-scores68 were previously calculated from ref. 22, in which each gene was calculated based on the frequency of observed nonsynonymous mutations (normalized to gene length and genome size) across 13 evolved lines23,24. G-scores for all genes in the genome of Dv and Mm were summed to get the “total observed G-statistic” and compared to the simulated “total expected G-statistic” by calculating a z-score as previously described22,68.

Scanning electron microscopy

Scanning electron microscopy (SEM) was used to visualize the overproduction of EPS from Dv-galUP32S using a 300 generation evolved line22,23,24. Monocultures were grown in 5 ml CCMA medium17 supplemented with 10 mM sodium acetate and 30 psi H2, whereas cocultures were grown in 20 ml CCMA medium supplemented with 30 mM sodium DL-lactate. Cocultures were incubated shaking at 37 °C until reaching early (OD600 = 0.11–0.15) and late stationary phase (OD600 = 0.4–0.45). 1 ml of culture was removed, diluted to an OD600 of 0.08 and fixed with glutaraldehyde (1% V/V) and sent for SEM to the Lawrence Berkeley National Laboratory.

Samples for SEM were prepared on 0.1% poly-l-lysine coated silicon wafers and a cell suspension was placed on the wafers for 15 min to allow cells to adhere. Excess cell suspension was then gently rinsed off with water and the wafers were fixed with 2.5% glutaraldehyde in Sodium Cacodylate (pH 7.2) for 1 h and post-fixed with 1% Osmium Tetroxide in Sodium Cacodylate (pH 7.2) on ice for another hour and finally rinsed three times with Sodium Cacodylate buffer (pH 7.2) to remove the remaining fixatives. Dehydration was performed through a graded ethanol series (20%, 40%, 60%, 80%, 90%, 100%, 100%, 100% at 7 min per step), followed by critical point drying using the Tousimis AutoSamdri 815 Critical Point Dryer (Tousimis, MD, USA). After the samples were dried, they were sputter coated with gold-palladium using a Tousimis Sputter Coater (Tousimis, MD, USA) to prevent charging in the microscope. Images were collected using the Hitachi S5000 Scanning Electron Microscope (Hitachi High Technologies America Inc, CA, USA).

Calculations of maximum growth rates and doubling times

To calculate the maximum growth rate of the planktonic and sediment phase we used the R package “growthrates”, which uses a logistic growth model, where μmax is the maximum growth rate and K is the carrying capacity69.

$$\frac{{dy}}{{dt}}={{{{\rm{\mu }}}}}_{\max }\times {{{\rm{y}}}}\left(1-\frac{{{{\rm{y}}}}}{{{{\rm{K}}}}}\right)$$
(1)

To calculate doubling time, we divided the ln(2) by the maximum growth rate.

$$\frac{{{\mathrm{ln}}}\left(2\right)}{{\mu }_{\max }}$$
(2)

Gene expression analysis

Paired-end Illumina reads were processed using TrimGalore version 0.4.363 and aligned to concatenated reference genomes of Desulfovibrio vulgaris Hildenborough (ASM19575v1 [https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_000195755.1/]) and Methanococcus maripaludis strain S2 (ASM1158v1 [https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000011585.1/]). Transcript abundance was quantified using kallisto v0.44.070 followed by separation of Dv and Mm reads prior to differential gene expression analysis using DESeq2 package v1.22.271 in R. Using kallisto transcript abundance estimates as TPMs (Transcript Per Million) for each organism, expression matrices were imported using the function “DESeqDataSetFromMatrix”. Specifically, differential expression analysis was performed by comparing planktonic expression to the sediment expression for all comparisons using the “collapseReplicates” function from DESeq2 for each day separately (e.g., Day 1 Planktonic vs Day 1 Sediment) and for all days combined (Days 1–6 Planktonic vs Days 1–6 Sediment) (Supplementary Data 3). Gene expression data of a wild-type (i.e., ancestral) Dv and Mm SynCom was obtained from a previous study21 (NCBI GEO database accession: GSE79022). The data were quantile normalized prior to performing principal component analysis to compare transcriptome states of the SynCom pre-evolution, after 1 K generations of laboratory evolution, and across the planktonic and attached phases over six days of culturing a minimal assemblage of the evolved community in the FBR.

Metabolic model refinement and integration

The genome-scale metabolic models of Dv (iJF744)46 and Mm (iMR539)47 were used in this study to represent Dv and Mm individually and further updated for accuracy. The published version of iJF744 metabolic model consisted of 1016 reactions and 951 metabolites. The gene-protein-reaction (GPR) relationship was represented using 744 genes of Dv in the network. These 1016 reactions were then curated as suggested in the publication. These updates included the renaming of the reaction ‘rxn05938’ as pyruvate synthase, removing the protein biosynthesis, DNA replication, and RNA transcription (i.e., reactions - ‘rxn13782’, ‘rxn13783’, and ‘rxn13784’) and fixing the capitalization of ‘Ex_cpd00218’ (Rxn name changed from ‘Ex_cpd00218(e)’ to ‘EX_cpd00218(e)’). Two reactions (‘rxn11934B_CC’ ‘QMO (Coculture)’ and ‘rxn11934B_SR - ‘QMO (Sulfate Reducing)’) were added in the place of one reaction (rxn11934) and also the exchange reactions of hydrogen, acetate and pyruvate were renamed to ‘EX_cpd11640(e)’, ‘EX_cpd00029(e)’, ‘EX_cpd00020(e)’ respectively. All of these updates were performed according to the ‘alterDvhModel.m’ script from the original publication. This curation of the iJF744 model resulted in 1014 reactions, 744 genes and 948 metabolites. We also updated the GPR of the reaction (rxn14403) as an “OR” relationship for the gene 208282 (DVU2776) because the sulfite reductase has been characterized to have multiple functions72,73. The published version of the iMR539 model consisted of 539 genes from Mm that were linked to a total of 668 reactions and 710 metabolites47. 51 exchange reactions in the Mm model were updated for the lower bounds to limit the metabolite exchanges for the syntrophic condition.

The integration of the two metabolic networks (iJF744 and iMR539) into a single syntrophy community model iSI1283 was achieved by interlinking the compartments of Dv and Mm through the known metabolite exchanges across Dv and Mm16. Initially, the metabolites of Mm were compartmentalized using standard tags for cytosol (“c0”) and extracellular (“e0”) compartments with ‘0’ at the end to specifically mention that these metabolites belonged to Mm. Hence, these tags identify Mm as “Organism 0”. Metabolite tags with (“c”) and (“e”) represent Dv. All the exchange reactions used for introducing metabolites to the extracellular compartment were standardized in “EX_{metabolite ID}[e0]” format. These exchanges were not removed because the ions and CO2 were allowed from extracellular compartment i.e., the media. The dependency of Mm on Dv was set by allowing interfacial constraints using the tDM [(t)ransport between (D)v and (M)m) reactions for hydrogen, formate and alanine that transports the metabolites from compartment [c] of Dv to compartment [c0] of Mm. Thus, the iSI1283 model integrated the two networks (iJF744 and iMR539) directly, given that the metabolites and reaction ids are the same but contains these additional tags to identify them.

After model integration, additional refinements to the individual Dv model were added to improve its proportionate predictions of hydrogen production from lactate oxidation by curating the model reactions from the iJF744 reconstruction46. The iSI1283 model was checked for elemental balance of key metabolites including reduced ferredoxin (cpd11620[c]), oxidized ferredoxin (cpd11621[c]), Menaquinol (cpd15499[c]), Menaquinone (cpd15500[c]), reduced DsrC (cpd18072[c]), oxidized DsrC (cpd18073[c]) that were initially proposed by ref. 46 in their iJF744 model for Quinone-interacting membrane-bound oxidoreductase (Qmo), Heterodisulfide reductase/NADH reductase (Hdr-Flx) and oligomeric lactate dehydrogenases (Ldh). According to ref. 46, for fermentative growth on a lactate medium, the overall expected stoichiometry is: Lactate + H2O -> Acetate + 2H2 + CO2. While checking for the elemental balance and expected stoichiometry, we have identified that the reduced ferredoxin (cpd11620[c]) and oxidized ferredoxin (cpd11621[c]) in one reaction (rxn06874) may have been misrepresented to carry 1 electron, hence we updated the reaction to represent the elemental balance to reflect these two metabolites as 2 electron carriers. This was performed using the COBRA toolbox script [model = addReaction(model,‘rxn06874’, ‘reactionFormula’,‘4 cpd11620[c] + 16 cpd00001[c] + 16 cpd00002[c] + cpd00528[c] <=> 16 cpd00009[c] + 2 cpd00013[c] + 4 cpd11621[c] + 6 cpd00067[c] + 16 cpd00008[c] + cpd11640[c]’,‘geneRule’,’(209712 and 209711 and 209708)’);]. In addition to this, reaction rxn14407 which accounts for periplasmic hydrogenases was also updated to reflect the representation of Cytochrome c2 (cpd00110[e]) and Cytochrome c3 (cpd00109[e]). This was also performed using the COBRA toolbox script [model = addReaction(model,‘rxn14407’, ‘reactionFormula’,‘2 cpd00110[e] + 2 cpd00067[e] <=> 2 cpd00109[e] + cpd11640[e]’,‘geneRule’,’((207234 and 207235) or (207390 and 207391) or (208022 and 208023))’);]. We also identified 50 bidirectional reactions that were constrained to be unidirectional by having an upper bound of 0 in iJF744 model (Supplementary Data 4). We reassigned them as bidirectional by changing the upper bound to be 1000. This allowed the model to release acetate into the extracellular environment, which is more accurate to what is observed during experimentation.

The final iSI1283 model was performance tested by comparing it with the core Dv and Mm mutualism model from ref. 17 by setting the same constraints for lactate and Dv biomass. iSI1283 predicted a ratio of hydrogen produced by Dv as 1.8-fold with respect to lactate; and methane release by Mm as 0.4-fold with respect to lactate; with an accuracy of 94% compared to the Stolyar et al., core model17. All the metabolic model files are available for download from our GitHub repository https://github.com/baliga-lab/Metabolic-Model-of-Syntrophy-for-Dv-and-Mm74 and are also provided in the Supplementary Data 4.

Flux predictions using state-specific models and analysis

The iSI1283 metabolic model was contextualized using gene expression data from planktonic and sediment phase obtained for day 1 to day 6 (generating 12 models). This step involved the use of the GIMME algorithm48. The GIMME algorithm takes the gene expression profile and reduces the model into a reduced network with all reactions that have genes expressed above the threshold. We applied the constraint-based method for simulating the metabolic steady-state of the SynCom model using flux-balance analysis (FBA)48,75. The initial validation steps involved checking the capacity of the syntrophic model to produce biomass in a defined media for the syntrophic coculture and checking whether Mm can produce methane in a Dv-dependent manner. In silico flux predictions were performed using the COBRA toolbox “optimizeCbModel” function and “fluxVariability” function in MATLAB. For predicting the day-wise phase-specific flux values, the 12 models (6 for planktonic and 6 for sediments) were further constrained by using Mm biomass by setting an upper bound for Mm that limits the Mm’s growth according to the phases and the day i.e., planktonic or sediment based on the experimentally measured Mm abundances. We used our models to predict key parameters of syntrophic growth, such as the biomass for two organisms (Dv biomass, Mm biomass), lactate exchange, acetate exchange, hydrogen production (used and unused), and methane production. Constraint-based modeling requires that one of the key parameters to be a measured value in order to predict the 5 other interplaying parameters. Here, we constrained the Mm biomass and then predicted the flux states by using either the Dv biomass or the methane as objective function. Final contextualization of the models integrated the measured ratio of Dv and Mm as weights to generate the most accurate predictions of flux states and microenvironment measurements of lactate, hydrogen and methane. All model simulations related to FBA were performed on the MATLAB_R2023b platform using the recent version of COBRA -The COnstraint-Based Reconstruction and Analysis toolbox (version 3.0)76. The GitHub repository https://github.com/baliga-lab/Metabolic-Model-of-Syntrophy-for-Dv-and-Mm74 contains all the details of the metabolic models. Supplementary Data 4 contains detailed spreadsheets related to iSI1283.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.