Synergistic epistasis enhances the co-operativity of mutualistic interspecies interactions

Early evolution of mutualism is characterized by big and predictable adaptive changes, including the specialization of interacting partners, such as through deleterious mutations in genes not required for metabolic cross-feeding. We sought to investigate whether these early mutations improve cooperativity by manifesting in synergistic epistasis between genomes of the mutually interacting species. Specifically, we have characterized evolutionary trajectories of syntrophic interactions of Desulfovibrio vulgaris (Dv) with Methanococcus maripaludis (Mm) by longitudinally monitoring mutations accumulated over 1000 generations of nine independently evolved communities with analysis of the genotypic structure of one community down to the single-cell level. We discovered extensive parallelism across communities despite considerable variance in their evolutionary trajectories and the perseverance within many evolution lines of a rare lineage of Dv that retained sulfate-respiration (SR+) capability, which is not required for metabolic cross-feeding. An in-depth investigation revealed that synergistic epistasis across pairings of Dv and Mm genotypes had enhanced cooperativity within SR− and SR+ assemblages, enabling their coexistence within the same community. Thus, our findings demonstrate that cooperativity of a mutualism can improve through synergistic epistasis between genomes of the interacting species, enabling the coexistence of mutualistic assemblages of generalists and their specialized variants.


Strains and Culture Conditions.
All the strains, culture conditions and the setup of the laboratory evolution experiment were the same as described before (Hillesland andStahl, 2010, Hillesland et al., 2014). Briefly, two clones of Desulfovibrio (Dv) and Methanococcus (Mm) were paired to setup 24 ancestral cultures in coculture medium A (CCMA) (Stolyar et al., 2007) under anaerobic conditions (80% N2:20% CO2 headspace) in Balch tubes. Cocultures were propagated weekly into a fresh media through 100fold dilutions and incubated either upright without shaking or in a horizontal position with constant shaking at 300 rpm. Laboratory evolution experiment was continued for 152 weeks and populations were archived as frozen glycerol stocks after generations 100, 300, 500, 780, and 1000 generations. Biomass collection was done as described before (Hillesland et al., 2014).

Sequencing of Evolved Cocultures.
DNA sequencing was performed for 13 of 22 evolved cocultures after 1000-generation and 9 evolved cocultures were sequenced at 100, 300, 500 and 780-generations. In addition, three End-Point-Dilutions from UE3 (EPD-03, EPD-09 and EPD-10) and HR2 (EPD-01, EPD-05 and EPD-10) cocultures evolved for 1000 generations and 3 clones of Dv and Mm from each of these EPDs were sequenced. For each sample, DNA was extracted with Epicentre Masterpure Kit (Epicentre Catalog number: MC85200). Sample and sequencing library preparation was done by using the Nextera DNA library preparation kit (Illumina) according to the manufacturer's instructions. DNA sequencing was performed in an Illumina Hiseq (generations 100, 300, 500, and 780) with 100 bp paired end sequencing or in an Illumina MiSeq sequencing instrument in the paired-end mode producing 2x250 bp long reads as described before (Hillesland et al., 2014).

Identification of Mutations in Evolved Cocultures.
Mutations accumulated in populations were determined by using a custom sequence alignment and variant calling pipeline (https://github.com/sturkarslan/evolution-of-syntrophy). This pipeline included quality control and trimming of the raw sequencing reads in fastq format by using Trim Galore software (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) and M. maripaludis (Genbank assembly: GCA_000011585.1) genomes and subsequent processing steps before calling the variants was done by following The Genome Analysis Toolkit (GATK) (DePristo et al., 2011) best practices. Briefly, reads were first aligned to the reference genome using Burrows-Wheeler Alignment Tool (bwa) (Li and Durbin, 2009) (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.9 (Li et al., 2009). 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

Single Cell Sequencing
For single cell sequencing, EPD-03 and EPD-09 from UE3 evolved cocultures were grown to midlog phase. Single cells of Desulfovibrio or Methanococcus were sorted into wells of a 96-well plate containing 3 µl of PBS and Buffer D2 from Repli-G single cell kit (Qiagen) by using Influx flow cytometer (BD). For each EPD, one plate for each of Desulfovibrio and Methanococcus was prepared. In order to lyse the cells, a freeze-and-thaw cycle was performed by first spinning the plates and freezing them at -20°C followed by thawing and re-spinning. Whole Genome Amplification (WGA) from single cells was performed by using REPLI-G Single Cell kit (Qiagen) according to manufacturer's instructions. We screened single amplified genomes (SAGs) with 16S universal primers for Desulfovibrio or Methanococcus to identify percentage of wells that did not contain any amplified product due to missing cells or failed WGA reaction. Wells with confirmed amplification were further treated with AmpPure XP magnetic beads (Beckman-Coulter) to clean and purify SAGs. A subset of SAGs was also analyzed with Bioanalyzer to confirm the size of the amplified fragments. Concentration of the SAGs passing the quality controls were determined by using Quant-iT PicoGreen dsDNA assay kit (Thermofisher). Nextera XT Library preparation kit (Illumina) was used for sample and library preparation for sequencing.
Sequencing was performed in HiSeq platform (Illumina) by using High-Output flow cell in 2x150 bp paired-end format. Sequence analysis including quality controls, trimming, alignment and variant calling was performed as described above.

Single Cell Lineage Tree Building
Variants identified from single cells of Desulfovibrio or Methanococcus for EPD-03 and EPD-09 of 1000-generation evolved UE3 cocultures were converted into a binary mutation matrix where each row was a unique mutation and each column was a single cell. Values in the binary mutation matrix were either 1 (mutation was seen in that particular cell), 0 (mutation was not observed) or 3 (there wasn't enough confident reads to assign the mutation). A variant was considered for the analysis only if its frequency was over 80% and was seen in at least two single cells. Mutation histories of single cells were determined by using SCITE algorithm (Jahn et al., 2016) with parameters -r 1 -l 900000 -fd 6.04e-5 -ad 0.21545 0.21545 -cc 1.299164e-05. SCITE used stochastic search to find the Maximum Likelihood tree of mutation histories in Newick format, which was converted to Cytoscape format for visualization purposes. This tree represents the predicted temporal order of the mutation events. Mutations were re-ordered by using information from the sequencing of the early generation cocultures if the order of the mutations couldn't be determined from the single cell mutational profiles due to noisy and missing data. Mutation tree was further annotated with gene functions, type of mutations and status of the mutations in early generations, and clonal isolates.

Calculation of G-scores
Based on the frequency of observed mutations (normalized to gene length and genome size) across 13 evolved lines, we calculated a G-score ("goodness-of-fit") to assess if the observed parallel evolution rate was higher than background as described before (Tenaillon et al., 2016).
Briefly, expected number of mutations ( ! ) for each gene in the genome was calculated as: where "#" is the total number of mutations, ! is the length of the gene and "#" is the total length of the coding genome. G-score for each gene ( ! ) was calculated as: where ! is the number of nonsynonymous mutations observed for gene across all evolved lines. G-scores for all genes in the genome of each organism were summed up to get the "total observed G-statistic" ( #%& ). In order to get the "total expected G-statistic" ( $'( ), we simulated Ntot number of mutations randomly across the protein-coding genome and calculated the mean and standard deviation of G-statistic from all simulations. We compared the observed and expected G-statistics by calculating a Z-score as follows; Where ( $'( ) and ( $'( ) are the mean and standard deviation of the G-statistics from 1000 simulations, respectively.

Density dilution assay
Ancestor and EPD cocultures were revived anaerobically in 18×150-mm balch tubes (Chemglass Life Sciences: CLS420901) from freezer stocks through dilution into CCMA media to ensure syntrophic growth and prevent carryover of glycerol into the fresh growth medium. EPD batch cultures (10 mL) were grown in anaerobic conditions with 80%:20% N2:CO2 headspace at 30°C without shaking. Growth of the cultures were monitored using a spectrophotometer (Spectronic 200: Fisher Scientific) to measure the optical density at 600 nm (OD600), measurements were typically taken twice a day until the cocultures reached stationary phase (~0.7 for EPD-03 and EPD-09). The cultures were kept in stationary-phase for approximately 20 hours in order to ensure similar growth phases before dilution into 96-well plates (Thomas Scientific: 1154Q44) for the density dilution assay. Stationary-phase Ancestor and EPD cultures were all diluted to the same starting optical density (as measured by the 96-well plate reader: BioTek EPOCH2T), so that all cultures in column 12 of each plate have the same starting optical density. Thereafter, each column between 11 and 2 received a volume of cells that would equate to a 1.5-fold dilution of the previous column and starting with column 11. The 1 st column of each 96-well plate contained only media as a control to identify potential contamination. All dilutions of cocultures into plates for the density dilution assay were done inside a Coy Anaerobic Chamber with an approximate atmospheric ratio of 95%:5% N2:H2. Plates were sealed with optically clear strong adhesive PCR films (115x100mm, Thomas Scientific: 4ti-0500/8), and the edges of these seals were coated twice in clear acrylic to further inhibit potential gas diffusion into the wells. Following inoculation, sealed plates were incubated at 30°C within the anaerobic chamber. Plates were removed from the anaerobic chamber twice per day to take growth measurements in the BioTek plate reader.
Plates in the plate reader were shaken linearly for 5 seconds prior to OD600 measurements at 30°C. During the transfer of plates from anaerobic chamber to plate reader, they were insulated between two 6-well plates filled with H2O that were also incubated at 30°C in order to maintain constant temperatures during plate transport from anaerobic chamber to plate reader. Density dilution assays were carried out for approximately ~4-5 days. A moving average with a window of two was applied to ODs for each density dilution assay to smooth the timeseries data. In order to establish a baseline for growth, a threshold was calculated for each EPD and ancestor strain based on the minimum carrying capacity from the first two dilutions (n=16). If cells from a well did not achieve an OD of that minimal threshold or higher, those cells were considered as not grown.

Clonal isolate pairings and measurement of Growth Rate and Yield.
Isolates of Desulfovibrio from EPD-03 or 09 of line UE3 were revived anaerobically using 10 ml of CCMA containing 10 mM:7.5 mM sodium lactate:sodium sulfite or 30 mM:20 mM sodium lactate:sodium sulfate, respectively, in balch tubes flushed with 80%:20% N2:CO2. Isolates of Methanococcus were revived anaerobically using 5 ml of CCMA containing 10 mM sodium acetate in balch tubes pressurized to 30 psig with 80%:20% H2:CO2. After the second transfer of revived isolates on their respective media, 0.1-0.2 ml of stationary phase cultures were combined in 20 ml of CCMA containing 30 mM sodium lactate in balch tubes flushed with 80%:20% N2:CO2.
After the second transfer, cocultures were stored as freezer stocks for future growth analysis.
Cocultures for growth analysis were revived anaerobically using 20 ml of CCMA containing 30 mM sodium lactate in balch tubes flushed with 80%:20% N2:CO2. Cultures were incubated at 37°C and shaken horizontally at 300 rpm. Optical densities (OD600nm) were monitored to assess growth and growth parameters were estimated using the fitting package grofit (Kahm et al., 2010).

Excess over Bliss analysis for measuring synergy
We adapted the Bliss Independence model (Borisy et al., 2003)  A clonal isolate pair combination for which EOB ≈ 0 has an additive behavior, whereas a pair with positive or negative EOB values has synergistic or antagonistic behavior, respectively. Error bars were computed by propagating the standard deviation of fractional effects.