The study of evolution in organisms with high mutation rates requires detailed data on the genomic composition of the population. Here we report on an innovative use of high-throughput sequencing technology combined with rapid experimental evolution that allowed us to follow the genetic diversification of four independent bacteriophage populations for 412 generations in extraordinary detail. One of our goals was to assess phylodynamic methods that use genetic information to infer key parameters of the population history. Because phylodynamic methods assume neutral evolution, we did not actively impose any selection pressures in our evolution experiment (in addition to those inherent in the bacterial cell environment). Tracing over 80′000 bacteriophage genomes, representing 884 distinct genotypes, we find that the patterns of viral diversification result in largely non-congruent genotype distributions, but also document multiple instances of parallel evolution. Unlike previous observations of parallelism in viral evolution, we do not just see the parallel evolution of mutations at single sites, but on the level of full-length genomes. Computer simulations that recapitulate our experiments in great detail show that this degree of parallelism is inconsistent with neutral evolution. We further show that the observed extent of parallel evolution biases phylodynamic analysis of migration rates, and wrongly estimates significant migration between the independent evolution lines. Our approach is applicable to many other viruses and paves the way to advanced population genetic investigations that aim to shed light on phenomena that relate to genetic linkage across the genome.

Genetic sequences of organisms contain imprints of their population biology and migration history. To recover this history, phylodynamic inference has been developed that combines phylogenetics with population biological models. Although phylodynamic methods are being applied widely, they have not yet been experimentally validated. To do this, we evolved a bacteriophage in four independent evolution lines for over 400 generations. We did not intentionally impose any selection pressure in these evolution experiments because phylodynamic inference works best in neutrally evolving systems. We found that often the same mutations and combinations of mutations evolve in the independent evolution lines. As a result, phylodynamics wrongly estimates migration between the lines. We conclude that the effects of selection should be included in phylodynamic inference.

Generating data on the evolutionary dynamics with long-term evolution experiments in a liquid handling platform

To generate data on the evolutionary dynamics, we evolved φX174 in 4 replicates in a constant and controlled environment of its host cells—Escherichia coli C (see Fig 1 and Materials and methods). During the evolution experiment, we did not intentionally impose any selection pressures, for example, by providing unusual host cells to which the phage was not adapted to, or by inhibiting phage replication. We could not rule out that selection pressures were operating though, and provide evidence for selection below. The rationale for not intentionally applying selection was to obtain insights into near neutral bacteriophage diversification and to generate a data set for phylodynamic analysis that conforms to the underlying assumptions of phylodynamics. We seeded each evolution line with the same, genetically heterogeneous stock population of φX174. We then serially transferred 2% of the bacteriophage population into wells containing fresh host cells for 412 generations. We transferred approximately every 26 min, consistent with the duration of one replication cycle of φX174. With a newly developed experimental protocol (see Materials and methods), we then extracted the DNA of 5% of the final bacteriophage population (generation 412), of an intermediate population (generation 196), as well as the original, “ancestor” bacteriophage population (generation 0), and sequenced this DNA using PacBio sequencing [15].

PacBio sequencing requires double-stranded DNA that is turned into a circular template. Because the average read length of PacBio is several times longer than the genome of φX174—over 38 kb in our experiment—the majority of reads contain multiple, (re-)sequenced copies of each sampled phage genome. By considering the consistency with which each nucleotide has been called at a given position, one can derive a quality score for each position, and generate the so-called circular consensus sequence [6,16]. Thus, the high error rate of PacBio, and other long-read sequencing platforms, of a single nucleotide read (approximately 10%–15%) is substantially reduced by reading each position multiple times [16]. To further increase the accuracy of our whole-genome sequences, we restrict the circular consensus sequences to quality scores higher than 47, which translates to a sequencing error rate of per nucleotide, or to less than one mis-called nucleotide per genome.

After these procedures to ensure sequence quality, we obtain over 110,000 reads covering the entire genome of φX174, representing over 31,000 distinct genotypes (see Fig 2). The majority of these genotypes are only present at very low frequency though. This raises the question whether they were present in the sample or whether they arose due to mutations introduced during the PCR or cell amplification involved in extracting the DNA and preparing it for sequencing, or due to errors introduced during PacBio sequencing that were not corrected through the circular consensus sequences.

Simulations guide thresholds to exclude sequences generated by DNA extraction or sequencing errors

In order to quantify the extent to which mutations are introduced during the DNA extraction and preparation, we performed in silico simulations. These simulations captured the population genetics of the serial passaging experiment, as well as the methodological steps involved in sampling, preparing and sequencing the viral populations. The exact steps of the experimental setup incorporated in the simulations are shown in Fig 3A and are explained in detail in the Materials and methods. In brief, the serial passaging experiment is simulated using a previously described sequence evolution simulation tool [17], which includes the replication and genetic diversification of the viral population as well as the bottlenecking involved in passaging and has been adapted to φX174. We also incorporated fitness effects of mutations in φX174 that have been experimentally estimated by Domingo-Calap and colleagues [18]. The simulations also accounted for the methodological steps involved in sampling, preparing and sequencing the viral populations. The sampling of the viral population and the DNA extraction are implemented as a series of bottlenecking and amplification steps that are involved in the protocol. Bottlenecks can lead to the loss of genotypes, while amplification steps can introduce mutations and distort the relative frequencies of genotypes in the sample. Lastly, the simulations incorporate the introduction of sequencing errors at a rate of per nucleotide (consistent with restricting the sequencing data to reads with a PHRED score of 47 or higher). The detailed simulation allows us to understand in which exact step sequence variation is introduced, and to which extent this contributes to the variation observed in the final reads.

The simulations allow us to estimate how frequently mutations introduced during the DNA extraction protocol generate genotypes that are not present in the evolved populations. We find that such genotypes with errors do not reach high frequency. Most of the genotypes with errors (70%–80%) are introduced during sequencing, while the different DNA extraction steps account for at most 30% of these genotypes. In none of our simulations did a genotype with errors occur more than three times (see Fig 3D).

Because our simulation showed that a genotype that occurs more than three times in the sequence data is unlikely to be the result of methodological error, we constrained our subsequent analysis to genotypes that are above this threshold. Of the 31,229 distinct genotypes in the sequencing data only 884 genotypes occur more than three times, representing 80,886 individual genomes. In all samples pooled over every evolution line at both time points and the ancestor population, there are 184 sites in the φX174 genomes that are mutated compared to the consensus. In the ancestral population, there are 53 genotypes and 38 mutated sites. All these statistics are calculated in S1 File. A list of all 184 mutations and their characteristics can be found in S1 Table.

Whole-population, whole-genome data show the evolution of mostly new genotypes in each evolution line

The sequencing data that we generated in our evolution experiments is so detailed—especially due to the fact that we have thousands of whole-genome sequences at each sampled time point—that it presents a challenge for visualization. To overcome this challenge, we present several, partly newly developed ways to visualize our whole-genome, whole-population sequencing data.

A conventional summary plot of the sequencing data for evolution line 1 at generation 412 is shown in Fig 2. This plot shows the frequency of non-consensus nucleotides at each position in the genome. In the sample shown, polymorphisms more frequent than 1:1,000 occur in 260 sites across the genome, of which the most frequent occur in the F gene that encodes for the major capsid protein (Fig 2A). This visualization of the sequencing data, however, ignores the linkage of mutations between different sites, and does not show the frequencies of each genotype.

Fig 2B displays the information on the diversity and linkage in a bacteriophage population sample by showing stacks of genotypes. This is similar to highlighter plots as used in the context of HIV, e.g., by Keele and colleagues [4]. The height of these stacks is proportional to the frequency of the genotype in the sample. The genotypes are displayed as a reduced genome that contains only the sites that are polymorphic. (In this regard, this plot differs from the highlighter plot that typically displays the entire genome.) This visualization, while containing information on the frequency of each genotype, does not reveal low-frequency genotypes. The sample of the virus population in evolution line 1 at generation 412 that we display in Fig 2B consists of 6187 genomes, representing 93 distinct genotypes. The most prevalent genotype was observed 2,035 times, the second most prevalent—the consensus in the ancestor population—was found 973 times. Twenty-three genotypes in this sample were found 4 times (the lowest frequency above the cut-off defined above). These low-frequency genotypes are not visible in Fig 2B.

To be able to visualize the full frequency distribution of genotypes in a sample, we developed what we refer to as genotype count distribution (Fig 2C). In this plot, distinct genotypes are indexed in order of decreasing frequency, and the logarithm of their frequency is shown as a histogram. This clearly displays the number and frequency of each genotype in a sample, but, unlike a highlighter plot, does not display the specific mutations each genotype carries.

Using these genotype count distributions, the diversification of the bacteriophage population over time can be visualized in great detail (Fig 4). The grey genotype count distributions at the left of each evolution line in Fig 4 show the ancestor population, which is used to seed every evolution line and consist of 53 distinct genotypes. The lightly colored genotype count distributions in each line show the composition of the bacteriophage populations at generation 196. They largely do not overlap with the ancestor population, indicating substantial shifts in the viral composition over time in each evolution line. There is also almost no overlap between the genotype count distributions between the evolution lines. Thus, most genotypes that evolve are new and unique to the evolution line in which they evolved—as would be expected in a regime of neutral evolution.

An alternative visualization of the evolution of the phage in the four lines is the genotype network shown in S3 Fig. While this graph shows only a subset of the 50 most prevalent genotypes, it clearly illustrates the independent evolutionary trajectories in each line.

Extensive parallel evolution of single mutations, combinations of mutations, and even whole-genome genotypes

Although we generate mostly unique genotypes over time in each evolution line in our passage experiment, there are genotypes that detected in multiple evolution lines (indicated by red circles in Fig 4). It is very unlikely that these genotypes have been unintentionally carried over into multiple evolution lines by cross-contamination during passaging the bacteriophage. This is because, first, on the 96-well plate we used for passaging, we separated the wells in which the bacteriophage was passaged by control wells that were seeded with bacteria to detect spill-over (see Materials and methods). Most of these control wells were not infected by spilling-over bacteriophage, and in case of infection, the experiment was restarted (see Materials and methods), making it unlikely that spill-over occured across two rows on the 96-well plates. Second, many of the genotypes that emerged in multiple lines were even further apart on the 96-well plate (quite a few between evolution lines 1 and 3, as well as 1 and 4 - see Fig 4). Lastly, the parallel genotypes are never the most frequent genotypes in any of the evolution lines in which they appeared (Fig 4). Thus, these genotypes most likely represent instances of parallel evolution.

It is important to emphasize that these instances of parallel evolution represent the parallel evolution on the level of whole genomes, rather than parallel evolution of single mutations. In total, we observe 26 genotypes that evolve in more than one evolution line. Compared with the ancestor genotype that is genetically closest, they acquired 1–5 additional mutations (mean 2.5). In total, these 26 genotypes carry 18 mutations that are not found in the ancestor population (see S1 File). Of these mutations, 33% are synonymous, compared to 19% in the ancestor, and 25% in the entire dataset. Establishing parallelism on the whole-genome level requires the whole-population, whole-genome sequencing approach presented here that also captures genotypes of low frequency.

To investigate to what extent the observed parallel evolution is consistent or inconsistent with an accumulation of unbiased and selectively neutral mutations, we performed randomization tests (see Materials and methods). Specifically, we compared the frequency of the observed parallel individual mutations and combinations of mutations with the frequency predicted under unbiased neutral evolution. (Fig 5).

On the level of individual mutations, we found 137 that emerged in only a single evolution line, while 47 emerged in multiple lines. The fraction of synonymous mutations does not differ significantly between parallel and non-parallel mutations (see S1 File). Seven mutations were shared in all four evolution lines (Fig 5A—blue bars). Randomization tests show that there is excess parallelism in mutations the are found in three or four evolution lines, i.e., that this level of parallel evolution is not consistent with unbiased, selectively neutral mutations (Fig 5A—gray bars). This extent of parallelism is similarly high as in previous evolution experiments using the same bacteriophage [12,13,19] or other viruses [20,21].

The linkage information in our data puts us into the unique position to investigate the parallelism of combinations of mutations across the entire genome (Fig 5B–5G). We observed the emergence of combinations of up to five mutations in multiple evolution lines (Fig 5E). Randomization tests show that the extent of parallelism is significantly higher than that predicted by neutral evolution. Randomization tests with an effectively shorter genome (see Materials and methods) show that the parallelism of combinations of mutations cannot be fully explained by the excess parallelism of single mutations. This strongly suggests that the selective advantage of the mutations in combination must on average exceed that of the individual mutations involved, i.e., the mutations in many of the parallel-evolved combinations are positively epistatic. More direct evidence for positive epistasis between mutations could be gained from more fine-grained data on the frequency of genotypes over time. This, however, will require to sample and sequence more frequently than every 200 generations.

As for the parallel genotypes, the parallel mutations and mutation combinations are very unlikely to have arisen by cross-contamination during passaging the bacteriophage. The most compelling piece of evidence against contamination and in favor of true parallel evolution is that the parallel mutations are found on mostly different genomic backgrounds in different evolution lines. First, we investigated how diverse the genotypes are that harbor the 47 single parallel mutations in each evolution line. We determined that each mutation is found on multiple genotypes but that the majority of these genotypes differ across evolution lines. On average, only 2.5% of the genotypes are the same across the evolution lines in which a parallel mutation is detected. If we weight the average by the frequency of each distinct genotype, the fraction is 3.5%—similarly low. To visualize the low overlap of these genotypes across evolution lines, we have complemented this analysis by genotype network graphs for each of the 47 mutations (see S4 Fig). We also performed this analysis of the genetic background for combinations of 2, 3, and 4 parallel mutations. Similar to single parallel mutations, only a minority of the genotypes that harbored the combination were found in multiple lines. The fractions are 6.1%, 18%, and 27% for combinations of 2, 3, or 4 mutations, respectively. If we weight the average by the number of genotypes, the fractions are 6.5%, 13%, and 20% for combinations of 2, 3, or 4 mutations, respectively.

Whether these combinations arose by successive point mutations or by recombination is difficult to ascertain. The multiplicity of infection (MOI) during the passaging experiment increased to 0.65, making it possible that bacterial cells were co-infected, enabling recombination. An analysis of genetic linkage (see Materials and methods), however, did not find any evidence for signatures of recombination.

Phylodynamic inference methods over-estimate migration rates between evolution lines

Viral sequences are widely used for phylodynamic inference which in turn has become key to study the spread and evolution of important viral pathogens between-host [22–26], as well as within-host [27]. However, parallel evolution can lead to biases in phylodynamic inference [28,29].

To investigate to what extent the parallel evolution we observed in our experiment biases parameters estimates, we conducted a phylodynamic analysis, estimating the migration rates between the independent evolution lines (see Materials and methods). Because the viral populations evolved independently in each line, there is no migration between the lines, and the estimated migration rates should not be significantly different from 0.

Fig 6A shows the distribution of estimated overall migration rates between evolution lines. In this distribution, the migration from the ancestor to the evolution lines has been subtracted. In a nutshell, the migration rate estimates are significantly different from zero, which they should not be because there is no migration between the evolution lines in the experiment. The point estimate is 0.00275 per generation, which corresponds to a per-generation probability of 0.275% for a genome to migrate from one evolution line to the other. Over the course of the 412 generation experiment, this migration rate amounts to a 68% probability for every single genome to migrate to a different evolution line. Qualitatively, these insights also hold when we exclude the 26 genotypes that evolved in parallel Fig 6B.

This positive migration rate, rather than reflecting real migration, constitutes a bias that arises because the analysis does not account for the selection pressures that generate the patterns of parallel evolution in our experiment. Thus, even in systems as ours, in which no intentional selection pressure is imposed on the viral population, selection is strong enough to bias phylodynamic inference.

To provide direct evidence for the role of parallel evolution in the over-estimation of migration rate, we conducted a phylodynamic analysis masking all mutations that appeared in multiple lines, i.e., parallel mutations, or just those that appeared in 3 or 4 lines which we call "excess parallel mutations" because only these are more numerous than neutral simulations suggest (see Fig 5A). We find that, in both case, the lower bound of the 95% confidence interval of the migration rate estimate is not above 0 anymore, meaning that the bias is not statistically significant anymore. This clearly shows that parallel mutations are the source of the overestimation bias.

Strains, growth conditions and bacteriophage titration

Except otherwise stated, Escherichia coli C (ATCC 13706) is cultured in LB Miller broth (Sigma-Aldrich, Cat. L3022) at 37 °C with aeration at 180 rpm. For infection with bacteriophage φX174 (ATCC 13706-B1, obtained from the DSMZ collection), the medium was supplemented with 10 mM MgCl₂ and 2 mM CaCl₂.

Amplification of φX174 was carried out by inoculating 50 mL of supplemented LB broth with 500 μL of stationary-phase E. coli and 250 μL of φX174 stock. Cultures were incubated overnight at 37 °C with aeration. Following incubation, the cultures were centrifuged to pellet bacterial cells, and the supernatant was filtered through a 0.22 μm filter to isolate the phage-containing lysate.

Bacteriophage titration on plates was carried out by plaque assay. Specifically, 6 mL of soft agar, LB Miller soft agar (LB broth containing 5 g/L agar (Sigma-Aldrich, Cat. A7002), 10 mM MgCl₂ and 2 mM CaCl₂) was mixed with 400 μL of stationary-phase bacteria and 20 μL of diluted φX174 lysate, then poured onto LB Miller hard agar plates from LB broth with agar (Sigma-Aldrich, Cat. L2897) and the same salt concentrations. Plates were incubated overnight at 37 °C in an inverted position, and plaques were counted to determine the initial phage concentration.

For liquid-based titration, a 96-well plate format was used. Each well was filled with 200 μL of LB Miller broth containing 1% (v/v) stationary-phase E. coli and φX174 diluted to an average of 0.5 plaque-forming units (PFU) per well. Plates were incubated overnight at 37 °C, and the presence of phage replication was assessed by spotting 2 μL from each well onto a fresh bacterial lawn prepared identically to those used in the plaque assay but without added phage. After a 5-hour incubation at 37 °C, wells with no evidence of phage replication were counted. The initial phage concentration was estimated by calculating the fraction of wells without detectable phage replication, assuming a Poisson distribution of phage particles across wells.

Experimental evolution of bacteriophage φX174

Bacteriophage phiX174 was experimentally evolved for 412 generations in 4 independent evolution lines. The experimental evolution of bacteriophage φX174 was carried out in 96-well plates using a Tecan Freedom Evo 200 automated liquid handling system (Tecan) with an integrated automated incubator (Liconic STX100, Liconic). At each transfer, 4 μl of supernatant was passaged into the next well.

The passaging terminated in the liquid handling system 3 times at generation 95, 150, and 309. At generation 95, only 2 μl of supernatant was passaged.

Bacteriophages were passaged from column to column on a 96-well plate approximately every 26 min, allowing for one complete replication cycle in each well [10]. On a single column, 4 replicated lines of bacteriophage φX174 evolve, each separated by a line containing only host bacteria. This allows the detection of contamination events: a delay in bacterial growth in these lines indicates contamination by the bacteriophage. For each generation n, we realise a transfer (described in the next paragraph) and then incubate the plates for 20 min at 37 °C and 95% relative humidity. Once a day, we make a glycerol stock by mixing 100 μL of 40% glycerol with 100 μL of the n-1 bacteriophage generation. The stocks are then stored at −80 °C. In case of bacteriophage contamination or other problems, the experiment is restarted using the last stock as a bacteriophage source.

Each transfer is composed of several steps with pipette sterilisation after each step (see below). First, on column n + 5, 5 μL of freshly grown overnight E. coli C are used to inoculate 200 μL of LB Miller supplemented with 10 mM MgCl₂ and 2 mM CaCl₂. Second, 4 μL of freshly mixed column n-1 is transferred into column n. As a result, bacteria are inoculated approximately two hours before bacteriophage infection: this assures that bacteriophage infection occurs on exponentially growing bacteria. In addition, this protocol aims to ensure that the MOI stays largely below one, which limits bacteriophage recombination. Third, if needed, glycerol stock of bacteriophages are made as previously described.

The entire experiment is carried out with reusable tips, which require sterilisation between each step to inactivate any bacteria or bacteriophage particle. To do so, the tips are washed twice in 5.6% bleach, flushed from behind with water, washed in 70% ethanol and then in milliQ water.

To determine the MOI, the densities of E. coli C and phage φX174 were quantified within the experimental wells. Host cell density was determined via standard colony-forming unit (CFU) assays. Briefly, 5 μL of an overnight E. coli C culture were inoculated into 200 μL of LB Miller medium per well. We incubated the bacterial populations at 37 °C for 120 and 135 min as the duration of the automated pipetting regime was slightly variable over the course of the experiment, and increased from 120 min at the start to approximately 135 min during the experiment. For each duration, samples were diluted and plated on Hard Agar plates in six replicates to back-calculate the precise density within each well. The bacterial densities were CFU/well and CFU/well for 120 and 135 min of incubation, respectively, corresponding to CFU/mL and CFU/mL.

Bacteriophage titers were quantified via PFU assays at generations 0, 196, and 412. Phage abundances are reported relative to a 4 μL volume, corresponding to the specific initiation and transfer volume propagated between wells each generation. Absolute viral abundances within this reference volume were PFU at generation 0, PFU at generation 196 and PFU at generation 412. Combining these PFU measurements with the number of bacteria per well, we obtain MOIs of 0.022, 0.31, and 0.65 at generation 0, 196, and 412, respectively.

Extraction of the φX174 DNA

Two protocols were developed to extract φX174 DNA from experimental samples, targeting different stages of the bacteriophage’s life cycle (dsDNA versus ssDNA) and using distinct amplification and extraction strategies. (New protocols were needed because PacBio long-read sequencing requires 1 μg of phage DNA, many orders of magnitude more than we have in each well during our experiment.) The dsDNA extraction protocol was applied to all samples, forming the basis of the main dataset.

To compare the protocols, both methods were applied to two samples. As shown in S2 Fig, the dsDNA protocol yielded a greater diversity and abundance of distinct genotypes. Due to its higher efficiency, it was selected as the standard extraction protocol for the study.

dsDNA extraction

A ready-to-use version of this protocol is publicly available on protocols.io (https://dx.doi.org/10.17504/protocols.io.q26g7q3e3lwz/v1).

φX174 DNA was extracted during the replication stage when the genome forms a double-stranded circular dsDNA intermediate. This replication intermediate was stabilized by chloramphenicol-induced inhibition of bacterial protein synthesis cycle [44,45].

First, high-titer lysates were generated by amplifying 10% of the phage population in each sampled well on agar plates. Plates were overlaid with soft agar containing 400 μL of stationary-phase E. coli and 20 μL of phage sample. After overnight incubation at 37 °C, phages were extracted by adding 3 mL of LB broth to the plate, gently mixing, and recovering 2 mL of the supernatant, which was then filtered through a 0.22 μm filter.

Second, 200 ml cultures of E. coli, infected with φX174 lysates, were incubated with chloramphenicol and harvested after 3 hours. To prevent recombination during phage replication, the MOI was kept below 1 (see validation below).

Third, DNA was extracted using a modified protocol based on the plasmid miniprep protocol (NEB, Cat. T1010), adapted for larger culture volumes and sequential spin-column loading to isolate phage-derived double-stranded DNA.

Last, the extracted DNA was then linearized using NEB Cutsmart Buffer (NEB; Cat. B7204) and PstI-HF restriction enzyme (NEB, Cat. R3140) and finally purified with AMPure XP (Beckman Coulter, Cat. A63880).

DNA quality and quantity were assessed using three methods: concentration was measured with the Qubit dsDNA BR Quantification Assay (Thermo Fisher Scientific, Cat. Q32850) on a Qubit Flex Fluorometer, chemical purity was evaluated with a Nanodrop spectrophotometer (Thermo Fisher Scientific), and fragment integrity and contamination were checked using an Agilent TapeStation with Genomic DNA ScreenTape (Agilent Technologies, Cat. 5067), following the respective manufacturer instructions.

To validate that the MOI remained below 1 during infection, phage and bacterial concentrations were monitored. Phage titers from two samples were determined after overnight amplification on agar plates and found to be PFU/mL and PFU/mL, respectively. Bacterial concentration after 1.5 hours of growth in the 200 mL cultures was measured at CFU/mL, yielding a total of CFU per flask. Given that no more than 2 mL of phage lysate was added per 200 ml flask, the resulting MOI was confirmed to be smaller than 1.

ssDNA extraction

A ready-to-use version of this protocol is publicly available on protocols.io (https://dx.doi.org/10.17504/protocols.io.261geqxmwg47/v1).

To isolate φX174 single-stranded DNA (ssDNA), we used a five-step protocol based on established precipitation, enzymatic digestion, and whole-genome amplification steps.

First, high-titer lysates were generated as described above.

Second, 1 mL of lysate was mixed with 500 μL of 30% PEG/1.5 M NaCl, incubated at 4 °C for 24 hours to precipitate the phage particles, and centrifuged twice to pellet the phages. The pellet was resuspended in DNase buffer (NEB, Cat. B0303S), treated with DNase I (NEB, Cat. M0303S) and RNase A (Thermo Fisher Scientific, Cat. EN0531) to remove host nucleic acids and heated to 75 °C for 15 min to inactivate the DNase. Capsids were lysed by incubating with 20% SDS and pronase (Sigma-Aldrich, Cat. P5147) at 40 °C for 2.5 hours.

Third, DNA was purified by ethanol precipitation. 0.5M NaCl and two volumes of cold ethanol were added, followed by overnight incubation at −20 °C. The DNA was pelleted, washed with 70% ethanol, dried, and resuspended in 100 μL of DNase-free water.

Fourth, whole-genome amplification was performed using Q5 polymerase (NEB, Cat. M0491) with φX174-specific primers (forward: 5′-ATC GCT TCC ATG ACG CAG AA-3′; reverse: 5′-TGC AGG TTG GAT ACG CCA AT-3′). Each 50 μL PCR reaction contained 21.5 μL nuclease-free water, 10 μL Q5 reaction buffer (NEB, Cat. M0491), 2 μL DMSO, 2.5 μL of each primer (10 μM), 1 μL dNTPs (10 mM), 0.5 μL Q5 polymerase, and 10 μL template DNA. Five parallel reactions were run per sample using the following program: 98 °C for 30 seconds; 32 cycles of 98 °C for 10 seconds, 68 °C for 30 seconds, and 72 °C for 108 seconds; followed by a final extension at 72 °C for 2 min. Successful amplification was confirmed using an Agilent TapeStation with Genomic DNA ScreenTape (Agilent Technologies, Cat. 5067).

Finally, PCR products were pooled, purified with 0.8 volume AMPure XP beads (Beckman Coulter, Cat. A63880), according to manufacturer instructions, and eluted in 80 μL of EB buffer (NEB, Cat. T1016) pre-warmed to 50 °C. DNA concentration and purity were assessed with a Qubit BR dsDNA assay (Thermo Fisher, Cat. Q32850) and a Nanodrop spectrophotometer.

PacBio sequencing

The DNA was sent for long-read sequencing to the Functional Genomics Center in Zürich. Sequencing was performed on a PacBio single-molecule sequencer Sequel SMRT Cell 1M.

The long-read sequence data generated in the context of this study is publicly available in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB80636 (https://www.ebi.ac.uk/ena/browser/view/PRJEB80636).

The supplementary zip-archive on zenodo (https://doi.org/10.5281/zenodo.22711242) contains a file “data/haplotypes_run2.csv” that contains the genotypes and their frequencies in all samples extracted from the long-read sequencing data. This dataset retains all the information in the aligned long-read sequences in an efficient way, and is the main dataset we used for our subsequent analyses.

Simulations of the passaging experiment, DNA extraction, and sequencing

The simulation of our entire experiment, entails first a simulation of the viral populations and its diversification throughout 400 passages. This results in approximately virions with potentially different genotypes (although many virions carry identical genomes).

Subsequently, the simulated virions are sampled binomially according to the volume sampled from the wells in our experiments, which leads to a loss of rare genotypes and a distortion of their original frequencies. The sampled virions are then plated, their DNA is extracted, and, when simulating the ssDNA extraction protocol (see above), PCR-amplified, or when simulating the dsDNA extraction protocol (see above), replicated in cells. The simulated plating, PCR amplification, and replication in cells can introduce errors into the genomes. Lastly, the simulated genomes are sequenced, which is implemented as binomial sampling and a worst-case sequencing error rate.

For all the details on these simulations, please consult the file “python-code/code/simulations.py” in the supplementary zip-archive on zenodo (https://doi.org/10.5281/zenodo.22711242).

Simulating sequences as expected under neutral evolution

To generate genotypes that would be expected under neutral evolution, we randomized the locations of mutations according to the following procedure.

First, we infer which of the mutations of a genotype have newly emerged during the passage experiment, rather than having been already present in the ancestral population. We do this by identifying the genotype in the ancestral population that has most mutations in common with the genotype in question. We assume this specific ancestral genotype to be the evolutionary origin, and consider all the mutations of the genotype not already present in the specific ancestor to have newly emerged. In case of multiple ancestral genotypes with the same number of shared mutations from the genotype in question, we assume the evolutionary origin to be the one that is most frequent in the ancestral population.

Second, we randomize the positions of these likely newly emerged mutations. Formally, if we have n newly emerged mutations, we sample without replacement n out of , where is the length of the genome of φX174. (Thus, if, for example, a genotype has likely acquired a mutation at positions 1,770 and 2,966, this randomization could generate positions 3,328 and 4,599.) For the sake of simplicity, we randomize only the positions, rather than the positions and specific nucleotides substitutions. As a result, parallel evolution becomes more likely in these simulations, and the resulting randomization test is conservative, i.e., it underestimates the extent and significance of a deviation from neutrality.

Third, we adopted this procedure to generate viral populations with mutational patterns that would be expected under neutral evolution for each evolution line. Repeating this 1,000 times for all 31,229 distinct genotypes, we can estimate the probability of a mutated position being shared across multiple evolution lines, and compare it to the observed frequency (Fig 5A). In addition to considering single mutated positions, we also investigated combinations of mutated positions (Fig 5B–5G). The clear conclusion from this analysis is that there is much more sharing of single mutated positions and combinations of mutated positions than expected under neutral evolution (see main text). This could be due to a selective advantage that these mutations confer, or to a higher mutation probability at the positions that are shared across the evolution lines.

Last, we investigated if combinations of mutations are shared more than expected given that already single mutations are shared more than expected under neutral evolution. To this end, we set the effective genome length to instead of . As a consequence, shared mutations become more likely. At this particular effective genome length, the observed number of shared mutated positions across four evolution lines is comparable to the neutral expectation (Fig 5A). The number of mutated positions shared across two or three evolution lines even exceeds this expectation, which makes the following analysis conservative.

In contrast to single mutated positions, the observed number of shared combinations of mutated positions is still much lower than the expectation at this low effective genome length (Fig 5B–5G). This is strong evidence that these combinations have a selective advantage that exceeds the potential advantage that the individual mutations confer, i.e., interact under positive epistasis. Note that, while we cannot disentangle potential mutational bias from selective advantage for single mutations, once the effect of single mutations is subsumed by shortening the effective genome length, deviations from the neutral expectation must be due to positive epistatic interactions between the mutated sites.

Linkage disequilibrium analysis to probe for genetic signatures of recombination

To investigate if recombination occurs to an extent that it leaves a detectable trace in the viral genomes, we performed a bioinformatic analysis of the sequencing data. Recent advanced methods to detect recombination are based on reconstructing phylogenies for multiple parts of the genome, for example GARD [46]. These, however, are challenging to apply because we have over 80,000 whole-genomes in our data set. For this reason, we specifically used older, more scalable methods that exploit patterns of genetic linkage that are indicative of recombination.

The well-known relationship between linkage disequilibrium and recombination rate by Hill and Robertson 1968 [47,48] cannot be used for our purpose because of some of the assumptions that go into the derivation of the relationship. The formula applies only if there are no new mutations (which occur during our experiment), if there is no selection (which we assume is acting), and, most importantly, when an evolutionary equilibrium has been reached (which is not the case in our experiments either).

We therefore used the even older formula for the decline of the linkage disequilibrium over time [49,50], similar to recent methods developed by Romero and colleagues [43] for the estimation of effective recombination rates of HIV within the infected individual. From one generation to the next, the linkage disequilibrium between two mutations, , declines by a factor . Hereby, p₁, p₂, and p₁₂ denote the frequencies of genotypes carrying mutation 1, 2, or both mutations, respectively, and c denotes the probability of recombination. The recombination rate is commonly assumed to scale with the distance between the mutations on the genome, and can hence be rewritten as c = c_1bp d, with c_1bp denoting the recombination rate at a genomic distance of a single base pair, and d denoting the genomic distance in base pairs.

Extrapolating this decline of D over time, we obtain the following relationship between the linkage disequilibrium D, the genomic distance of two mutations d, the recombination rate c_1bp, and the generation g in our experiments:

(1)

We fitted this relationship (Eq 1) to all the mutation pairs that have a non-zero D in the ancestor population. We normalized all measurements of D by dividing them by the non-zero D in the ancestor population. We obtain an estimate for c_1bp that is not significantly different from 0 (). Thus, there is no signal for linkage disequilibrium decline over generations, and hence no evidence for recombination in our experiments. While this statistical test does not rule out any instance of recombination, it supports the assumption that recombination has not affected the evolutionary dynamics in a major way.

Phylodynamic estimation of migration rates between the four evolution lines

To investigate to what extent the instances of parallel evolution that occur in our evolution experiment can bias the estimation of migration rates between the evolution lines, we conducted a phylodynamic analysis. As the evolution lines in our experiments evolve independently, i.e., there is no migration between these lines, the migration rates should be estimated to be zero.

In our experiments, we obtained a total of 80,886 whole-genome sequences that appeared more than 3 times in one of the samples. Of these sequences only 884 are distinct, i.e., many of the sequences are identical. This leads to the problem of polytomies in the reconstruction of the phylogenies, on which a phylodynamic analysis is built.

The high number of sequences and the polytomies make a one-step Bayesian estimation of the migration rates using the software BEAST [51–54] prohibitive. We therefore resorted to a simpler, non-Bayesian two-step estimation procedure, ignoring the phylogenetic uncertainty inherent in the reconstructure of the tree. Technically, we conducted this analysis using treetime [31], an efficient algorithm to generate timed phylogenetic trees from sequence data. Treetime allows a two-step mugration analysis [30], which estimates migration rates from a phylogenetic tree analogously to mutation rates by adding a locus to the genomes that has its own mutation rate and designates the deme in which the sequences was sampled. The mutation rates between the, in our case, five demes then represent estimates of the migration rate between them.

To circumvent having to analyze all 80′886 sequences, we adopt the common strategy to subsample. We subsample 900 sequences, 100 from the ancestor population and another 100 from generation 196 and 412 of each of the four evolutiom lines, and repeat this analysis 100 times. This subsampling approach captures parts of the phylogenetic uncertainty that we did not incorporate directly.

Another complication that arises in the phylodynamic analysis of our experimental data is that the ancestor population, with which the four evolution lines were seeded, consists of multiple genotypes that arose through a common evolutionary history. This “panmictic past” will leave a signal for migration even in the independent lines. Our solution of this issue is to assign the ancestor its own deme, from which migration to each evolution line is possible. (Each evolution lines is also assigned a deme of its own.) This way, the common evolutionary history of the genotypes in the ancestor population can be separated from the potential migration of genotypes between evolution lines in the analysis (see S2 Table).

S1 Fig. Population dynamic simulation showing that phage-resistant bacteria cannot evolve in our experimental setup.

If bacteria evolve that are completely resistant to phage and carry no resistance cost, they are washed within a after a few transfers. This even applies to levels of resistant bacteria that are unrealistically high, such as the hypothetical level of 30% used in the simulation shown: they are washed out within five transfers. This is because in our experiment we transfer 2% of the phages (and some bacteria) every 26 min, (which corresponds to their lysis time) and transfer them into a well that has been seeded with a billion bacteria from a stock. The Python code of this simulation can be found in the file “python-code/code/simulation-showing-the-impossibility-of-bacterial-phage-resistance.ipynb” in the supplementary zip archive on zenodo (https://doi.org/10.5281/zenodo.22711242).

https://doi.org/10.1371/journal.pbio.3004064.s001

(PDF)

S2 Fig. Comparison of the dsDNA extraction protocol and the alternative ssDNA extraction protocol in 2 of our 9 samples.

(A) Comparison of the frequencies of the genotypes. (B) Comparison of the genotype counts distributions. The top and bottom rows display the genotype counts order by decreasing counts according to the dsDNA and ssDNA extraction protocol, respectively. The dsDNA extraction protocol was more efficient as it generates more distinct genotypes (see S1 File). The data and code needed to reproduce this figure can be found in the supplementary zip-archive on zenodo (https://doi.org/10.5281/zenodo.22711242).

https://doi.org/10.1371/journal.pbio.3004064.s002

(TIFF)

S3 Fig. Genotype network plot of the 30 most frequent genotypes.

Each genotype is represented by a pie chart, with the size of the pie chart representing the total abundance of that genotype across all samples. The consensus genotype is labelled. The pie chart then displays how the total abundance is distributed between the different samples. The pie charts are connected if and only if their genotypes differ by exactly one mutation (single-nucleotide substitution, or an insertion or deletion). In general, the different evolutionary lines cluster together as they contain closely related genotypes. In particular, all the non-ancestral genotypes shown that appeared in evolutionary line 3 are at least 2 mutations different to all other non-ancestral genotypes shown. The data and code needed to reproduce this figure can be found in the supplementary zip-archive on zenodo (https://doi.org/10.5281/zenodo.22711242).

https://doi.org/10.1371/journal.pbio.3004064.s003

(PDF)

S4 Fig. Genotype network plots of genotypes containing specific parallel mutations.

For each parallel mutation (those that occur in at least two different lines), a genotype network plot is shown for all genotypes present during the experiment containing that parallel mutation. The networks are presented in the same way as in S3 Fig. There is a large variance in the number of genotypes containing a specific parallel mutation, though in general genotypes are found mostly in a single evolutionary line. This is evidence against contamination, and supports that these mutations are the result of true parallel evolution. The data and code needed to reproduce this figure can be found in the supplementary zip-archive on zenodo (https://doi.org/10.5281/zenodo.22711242).

https://doi.org/10.1371/journal.pbio.3004064.s004

(PDF)

S2 Table. Migration rate estimate from each deme to each other deme (average over estimates for the 100 subsamples).

Deme A is the ancestor population, demes EL1–4 are evolution lines 1–4. The direction of migration goes from column to row. The estimated rates of migration from the ancestor population to the evolution lines can be seen in the first column of the matrix, and have consistently high values. This is consistent with the fact that each evolution line has been seeded with the ancestor population, i.e., migration actually happened in the experiment. The other four columns show estimates of migration between the evolution lines, and from the evolution lines back to the ancestor. These migration rates are zero in the experiment, but are wrongly estimated to be non-zero. This bias arises because of parallel evolution between the evolution lines. Between evolution lines 1 and 2, for example, there are particularly many genotypes that evolved in parallel by generation 412 (see Fig 4), leading to the largest migration rate estimates between these lines (0.00168 and 0.00182). The basis to calculate the overall migration rate estimate between the lines were columns 2–5 in the migration rate estimate matrices. Leaving out the first column in calculating the overall migration rate ensures that the migration rate is not confounded by the “panmictic past” of the ancestor population. The data and code needed to reproduce this table can be found in the supplementary zip-archive on zenodo (https://doi.org/10.5281/zenodo.22711242).

https://doi.org/10.1371/journal.pbio.3004064.s007

(PDF)