Origin of biogeographically distinct ecotypes during laboratory evolution

Abstract

Resource partitioning is central to the incredible productivity of microbial communities, including gigatons in annual methane emissions through syntrophic interactions. Previous work revealed how a sulfate reducer (Desulfovibrio vulgaris, Dv) and a methanogen (Methanococcus maripaludis, Mm) underwent evolutionary diversification in a planktonic context, improving stability, cooperativity, and productivity within 300-1000 generations. Here, we show that mutations in just 15 Dv and 7 Mm genes within a minimal assemblage of this evolved community gave rise to co-existing ecotypes that were spatially enriched within a few days of culturing in a fluidized bed reactor. The spatially segregated communities partitioned resources in the simulated subsurface environment, with greater lactate utilization by attached Dv but partial utilization of resulting H₂ by low affinity hydrogenases of Mm in the same phase. The unutilized H₂ was scavenged by high affinity hydrogenases of planktonic Mm, producing copious amounts of methane. Our findings show how a few mutations can drive resource partitioning amongst niche-differentiated ecotypes, whose interplay synergistically improves productivity of the entire mutualistic community.

Fig. 1. Partitioning of attached and planktonic syntrophic communities in FBRs.

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.

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

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”. 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 (Mmp1718_p.R23fs) compared to the attached phase (Mmp1718_wt 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 galU_P32S (DVU1283); dashed line indicates linear regression and R² 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 galU_P32S mutation, but not by ancestral Dv-galU_wt cells, which was confirmed across three independent evolved lines. Source data are provided as a Source Data file.

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

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 log₁₀ 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.

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

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) H₂ production by Dv (Total bars: hashed and non-hashed), H₂ utilization by Mm (non-hashed), and (H) amount of methane produced by Mm. I A heuristic model for the predicted differences in relative flux (mMol⋅gDCW⁻¹⋅h⁻¹) between planktonic and attached cells and ecotypes (Dv-galU_wt —Dark orange: Mmp1718_p.R23fs —Light grey and no archaellum; Dv-galU_P32S —Light orange and thick membrane: Mmp1718_wt —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.