Synergistic epistasis enhances the co-operativity of mutualistic interspecies interactions

Distribution, frequency, and functional implications of mutations during laboratory evolution of obligate syntrophy

We evaluated whether the selection of mutations in the same genes (i.e., “parallel evolution” [17]) had contributed to improvements in syntrophic growth of Dv and Mm across independent evolution lines, all of which started with the same ancestral clone of each organism. The goal of this analysis was to focus on generalized strategies for adaptation to syntrophy, irrespective of the culturing condition so we investigated parallelism across both U and H lines. Based on the number of mutations (normalized to gene length and genome size) in Dv and Mm across 13 evolved lines (six lines designated U for “uniform” conditions with continuous shaking and seven H lines for “heterogenous” conditions without shaking), we calculated a G-score [18] (“goodness-of-fit”, see “Methods” section [18]) to assess if the observed parallel evolution rate was higher than expected by chance. The “observed G-score” was calculated as the sum of G-scores for all genes in the genome of each organism; mean and standard deviation of “expected G-scores” were calculated through 1000 simulations of randomizing locations of observed numbers of mutations across the genome of each organism. The observed total G-score for Dv (1092.617) and Mm (805.02) was significantly larger than the expected mean G-score (Dv: 798.19 ± 14.99, Z = 19.63 and Mm: 564.83 ± 15.95, Z = 15.06), demonstrating significant parallel evolution across lines.

With the exception of five high G-score genes (DVU0597, DVU1862, DVU0436, DVU0013, and DVU2394), which were mutated during long term salt adaptation of Dv [19], mutations in other high G-score genes appeared to be putatively specific to syntrophic interactions. Altogether, 24 genes in Dv and 16 genes in Mm associated with core processes had accumulated function modulating mutations across at least 2 or more evolution lines (Fig. 2 and Supplementary Table S1). Signal transduction and regulatory gene mutations (seven in Dv and six in Mm) represented 19.9% and 27.2% of all mutations in Dv and Mm, respectively, similar to long term laboratory evolution of E. coli [18], potentially because their influence on the functions of many genes [20, 21]. We also observed missense and nonsense mutations in outer membrane and transport functions (four genes in Dv and three genes in Mm). For example, the highest G-score gene in Dv, DVU0799—an abundant outer membrane porin for the uptake of sulfate and other solutes in low-sulfate conditions [22], was mutated early across all lines, with at least two missense mutations in UE3 (S223Y) and UA3 (T242P). Notably, the regulator of the archaellum operon (MMP1718) had the highest G-score with frameshift (11 lines) and nonsynonymous coding (2 lines) mutations [23]. Similarly, two motility-associated genes of Dv (DVU1862 and DVU3227) also accumulated frameshift, nonsense and nonsynonymous mutations across 4 H and 3 U lines. Together, these observations were consistent with other laboratory evolution experiments performed in liquid media [24], suggesting that retaining motility has a fitness cost during syntrophy [25, 26].

Consistent with our previous observation that obligate mutual interdependence drove the erosion of metabolic independence of Dv [5, 27], mutations in SR genes were among the top contributors to the total G-score in Dv (DVU2776 (74.7), DVU1295 (46.5), DVU0846 (42.9), and DVU0847 (22.3)). However, it was intriguing that DVU2776 (DsrC), which catalyzes the conversion of sulfite to sulfide, the final step in SR, accumulated function modulating but not loss-of-function mutations across six lines. The functional impact of these mutations is not clear but it is possible that these changes might alter previously suggested alternative roles for this gene, including electron confurcation for the oxidation of lactate [28], sulfite reduction, 2-thiouridine biosynthesis and possibly gene regulation [29].

Analysis of temporal appearance and combinations of mutations across evolution lines

Growth characteristics of all evolution lines improved by the 300th generation [4], and in some lines even before the appearance of SR− mutations, indicating that mutations in other genes had also contributed to improvements in syntrophy. Each evolution line had at least 8 and up to 13 out of 24 high G-score mutations in Dv, while Mm had mutations in at least 5 and up to 10 out of 16 high G-score genes. We interrogated the temporal order in which high G-score mutations were selected and the combinations in which they co-existed in each evolution line to uncover evidence for epistatic interactions in improving obligate syntrophy. Indeed, missense mutations in DsrC (DVU2776) were fixed simultaneously with the appearance of loss of function mutations in one of two sigma 54 type regulators (DVU2894, DVU2394) in lines HA2, and UR1 (P = 5.40 × 10⁻⁵). In rare instances, we also observed that some high G-score mutations co-occurred across evolution lines, e.g., two U- and one H-line consistently showed for at least two time points a mutation in DVU1283 (GalU) coexisting with mutations in DVU2394 (P = 5.04 × 10⁻³). More commonly, the combinations of high G-score gene mutations varied across multiple lines. In fact, no two lines possessed identical combination of high G-score gene mutations (Fig. 3A, B), and many high-frequency mutations were uniquely present or absent in different lines (Fig. 3C, D).

Mutations in high G-score genes appeared consistently in all evolution lines (P < 1.82 × 10⁻²), although in a different temporal order in each line. Mutations in the same high G-score genes appeared at different times (e.g., whereas mutations in SR gene DVU0847 first appeared in the 300th generation of HA2, they appeared much later in HR2 and HA3) (Fig. 3A). Similar temporal patterns of appearance and co-occurrence of mutations were observed in Mm (Fig. 3B). We also discovered evidence for temporally nested fixations, wherein prior to fixation of a mutation from an earlier generation, another mutation selected in a later generation gradually increased in frequency towards fixation e.g., DVU0799, DVU0001, and DVU1283 in HA3 (P = 1.19 × 10⁻³); and MMP1718 and MMP0335 in UA3 (P = 1.81 × 10⁻⁴). Moreover, there were many cases of simultaneous fixation of mutations in multiple genes (e.g., DVU0799 and DVU1214 in HE3; MMP0378, MMP0705, MMP0986, and MMP1170 in HA2) suggesting that hitchhiking may be common [30, 31]. However, given that samples were only sequenced every 250 generations, we cannot rule out the possibility of each mutation sweeping separately during that time interval. These observations lead us to conclude that mutations that were commonly selected may simply have additive effects on fitness, arising at different times in different populations because of chance (i.e., they became available for selection at different times in different populations).

The longitudinal analysis revealed a cross-species selection event that resulted in the replacement of the dominant clones of both species with new clones containing different mutations. Between generations 500 and 780 of HS3, the dominant Dv (harboring high G-score mutations DVU0799, DVU0597, DVU0797) and Mm (harboring dominant mutations in MMP1255, MMP1611, MMP1362, and MMP1511) clones disappeared (Fig. 3E). At the same time a new Dv clone with mutations in DVU2394, DVU2451, DVU1862, and intergenic region IG_184033, and a new Mm clone with mutations in MMP0952, MMP1077, and MMP1479 were selected (Fig. 3E). One explanation for this phenomenon is that, coincidentally, rare clones in both Dv and Mm acquired beneficial mutations and outcompeted dominant clones in the same 250 generation interval of evolution. Another possibility is that selection of a new dominant clone in one species changed the selection environment for the other, allowing its rare clone to take over. Which species might have started this process is unclear because there are no samples available in the 250 generations during which the sweep occurred. However, information about the new mutations, their functions, and parallel evolution may provide hypotheses. The novel mutations in DVU2394 (a sigma 54-dependent transcriptional regulator) and DVU2451 (a lactate permease) co-occurred in at least two lines including HS3 (P = 3.93 × 10⁻²) and appeared individually in only three other lines, suggesting that the two genes might be beneficial and also functionally coupled in the context of promoting syntrophy. Notably, we demonstrate later through single-cell analysis that mutations in DVU2394 occurred subsequent to mutations in DVU2451, but exclusively in the SR− lineage within UE3. Interestingly SR− mutations were never selected through 1000 generations in HS3 line. Conversely, while mutations in DVU2394, DVU2451, and DVU1862, all co-occurred in UE3, they did not sweep through the population, underscoring how improvements to syntrophic interactions occurred through multiple distinct trajectories in terms of the order and combinations of mutation selection. In other words, this cross-species selective sweep occurred only in HS3, suggesting one of several features unique to this line was responsible, including simultaneous selection of mutations in DVU2394, DVU2451, and DVU1862, the overall mutational landscape of HS3 between generations 500 and 780, or mutations unique to HS3. Interestingly, fixed mutations that were observed only in HS3 were in Mm (MMP0952, MMP1077, MMP1479) and their appearance coincided with the selective sweep between 500 and 780 generations. Regardless of the mechanism, it is interesting that a new mutation(s) in Mm appears to have selectively swept high G-score mutations across both interacting organisms, strongly suggesting that the new Mm genotype conferred a fitness advantage to a specific lineage of genotypes in Dv that were in low abundance prior to the sweep.

Characterization of evolutionary lineages and interspecies interactions in minimal assemblages at single-cell resolution

We performed EPDs from the 1 K generation of HR2 and UE3 lines to generate simplified sub-communities that represented minimal sets of genotypes with growth phenotypes comparable to the 1 K culture (see “Methods” section). While two EPDs from each line represented the dominant SR− subpopulation of the 1 K evolved line, we also recovered an SR+ subpopulation that co-existed within each line albeit at much lower abundance and below the detection limit of bulk mutation analysis of the parental culture (Fig. 4). Finally, we isolated evolved clones of each organism by streaking EPDs on agar plates containing nalidixic acid and neomycin, taking advantage of chromosomally integrated selection markers in Dv and Mm, respectively. Altogether, three clones of Dv and Mm from each EPD were isolated and re-sequenced. The distribution of unique mutations across EPDs and clonal isolates added evidence for coexistence of distinct lineages of one or both organisms within each evolved line. Logically, all high-frequency mutations in an asexual population must be linked on the same genetic background. As expected, all 15 high-frequency mutations detected in the 1 K generation of UE3 were present only in sub-communities with the SR− mutations (EPD-03 and EPD-10). By contrast, at least 11 mutated loci (ten genic and one intergenic) in the SR+ sub-community (EPD-09) were not detected in the SR− sub-communities or in 1 K bulk sequencing of UE3, demonstrating that the EPD-09 assemblage was made up of rare Dv lineages (Fig. 4A and Supplementary Table S2). Strikingly, both Dv and Mm lineages in the SR+ assemblage of HR2 were distinct from lineages in the SR− EPDs, and below detection limit in 1 K bulk sequencing (Fig. 4B, Supplementary Fig 1 and Supplementary Table S2). Thus, the existence of genotypically distinct subpopulations of Dv and Mm with the parental growth phenotype suggested that specific interactions across multiple evolved genotypes of the two organisms could have emerged during their syntrophic evolution.

We further investigated the evidence for specific interactions among evolved genotypes using single-cell sequencing of SR− (EPD-03) and SR+ (EPD-09) assemblages from UE3. We sorted, amplified, and re-sequenced the genomes of single cells of Dv (94 from EPD-03, and 94 from EPD-09) and Mm (87 from EPD-03, and 72 from EPD-09) to reconstruct lineages of both organisms within each EPD ([32] and “Methods” section, Supplementary Table S3). Using stringent cutoffs (fold coverage ≥ 8, number of cells with mutation ≥2, frequency ≥ 80%) and consensus mutation calling using varscan [33], GATK [34], and Samtools [35]_, we identified across single cells of Dv 16 of 17 and 3 of 12 mutations detected in bulk sequencing of EPD-03 and EPD-09, respectively. Similarly, we identified across Mm single cells seven of seven and six of seven mutations from bulk EPD-03 and EPD-09, respectively.

Using a mutation lineage inference algorithm SCITE [36] and cross-referencing with longitudinal sequencing data from 100, 300, 500, 780, and 1000 generations, bulk sequencing of EPDs, single-cell sequencing, and sequencing of clonal isolates, we reconstructed the lineage and timeline of mutations that shaped the evolution of syntrophy in SR− and SR+ communities within UE3 (see “Methods” section) (Fig. 5 and Supplementary Figs. 2–3). As expected, the two EPDs shared a core lineage of events that included sequential accumulation of high G-score mutations in the early stages of evolution in both organisms. While the Mm lineages across EPDs had few differences, lineages of Dv were strikingly different across the SR− and SR+ communities. The SR− mutations in the EPD-03 lineage were followed by selection of mutations in at least six regulators, and complex radiating branches with many coexisting sub-clones, suggesting that loss of SR in the EPD-03 line might have promoted the selection of mutations in regulatory genes. Altogether, the observation that dominant lineages were excluded in the minimal community assemblages of EPD-09, demonstrates the coexistence of distinct high abundance (SR−) and low abundance (SR+) lineages within the same evolved population (Figs. 4 and 5 and Supplementary Figs. 1–3). A surprising observation is that the SR+ clone that remained in the population subsequent to the evolution of SR− was not simply the dominant clone without the SR− mutation. Instead, it was a rare genotype with different mutations from the dominant population.

Investigation of cooperativity and synergistic interspecies interactions

We performed a density dilution assay to investigate the evidence for improved cooperativity due to interactions among specific genotypes of Dv and Mm within microbial community assemblages of the two EPDs [37, 38]. Briefly, we generated a dilution series of both EPD and ancestral cell lines in 96-well plates and determined growth rate, carrying capacity, and a minimal cell density that supported syntrophic population growth (See “Methods” section, Supplementary Fig. 4). Both EPDs could initiate growth at significantly lower cell density relative to the ancestral coculture. EPD-03 initiated growth at a 1.5-fold lower cell density with faster growth rate and lower carrying capacity relative to EPD-09, explaining how the two EPDs co-existed in vastly different proportions in UE3 (>80% EPD-03 vs, <1% EPD-09, Fig. 6A). These data make a compelling case for the emergence of increased cooperativity among Dv and Mm lineages during the evolution of syntrophy.

We investigated whether the increased cooperativity could have emerged through synergistic interactions between Dv and Mm lineages by characterizing individual and combined contributions of the two evolved partners towards improved growth characteristics. By comparing pairs of evolved and ancestral clones (Dv_EV × Mm_EV, Dv_Ac × Mm_Ev and Dv_Ev × Mm_Ac, Dv_Ac × Mm_Ac), we determined that each evolved clonal isolate had contributed individually to significant improvement in growth rate and yield (Fig. 6B and Supplementary Table S4). The improvements were maximal, and comparable to growth characteristics of the parental EPD, when both partners in the interacting pair were evolved clonal isolates (Dv_Ev × Mm_Ev). This result demonstrated unequivocally that increased cooperativity had emerged from synergistic interactions between the evolutionary changes in both species within each EPD, with proportional antagonistic effect on growth yield [39] (Fig. 6C). The higher growth rate of EPD-03 and higher carrying capacity of EPD-09 (both relative to the other EPD) gives mechanistic insight into coexistence of SR− and SR+ sub-communities as a r– and K-strategists, respectively (Fig. 6D, Supplementary Fig. 4). Notably, the few mutations that differentiate genotypes of each clonal isolate appear to manifest in variation in growth rate and yield, demonstrating that productivity of Dv_EV × Mm_EV interactions are genotype-specific, even within the same EPD (Fig. 6D, Supplementary Fig. 5, Supplementary Table S2).

Source: Ecology - nature.com