Love the one you’re with: replicate viral adaptations converge on the same phenotypic change

Adaptation experiment

ID11 (GenBank accession number AY751298; Rokyta et al., 2005) is a single-stranded DNA bacteriophage of the family Microviridae. It has a genome of 5,577 bases encoding 11 genes arranged in the same way as G4 (from which it differs by ∼3%). We used the nine first-step beneficial mutations obtained by Rokyta et al. (2005) via flask-passaging and adapted each of these genetic backgrounds in eight-fold replicate for 72 total lineages in multiwell plates for 100 passages. Passaging began by adding 10⁴ phage into 500 µL of E. coli C host cells at a density of 10⁸/mL. Passaging was done by allowing phage/cell growth for 25 min in flat-bottomed 48-well plates at 37°C in an incubator shaking at 200 rpm. Plates were then put on ice and 5 µL of each well’s volume was transferred to the wells of a fresh plate containing naive hosts in the same volume and density. During growth, wells were sealed with a double-layer of gas permeable membrane to reduce the likelihood of contamination between wells. Five transfers occurred each day and plate were stored overnight in the refrigerator.

Sampling and sequencing

After 100 transfers we plated and picked five isolates from each lineage. Whole genomes were then sequenced using the Fluidigm Access Array Platform (South San Francisco, CA) and Roche 454 Genome Sequencing FLX (454 Life Sciences, Branford, CT) using a novel tagging method so that many small genomes could be sequenced at once while retaining linkage information. Raw Roche 454 unclipped DNA sequence reads were cleaned, assigned to barcode and amplicon, and mapped to the ancestral genome using the R package rSFFreader (R Core Team, 2014; Settles et al., 2011) and a custom R script (See Supplemental Information for more detail). Mutations were accepted as real only when coverage at the site where a mutation was detected was ≥10 and the mutation was present in >90% of the reads (mean coverage, 46×). To avoid strong sample size effects, we removed wells from the analysis with three or fewer successfully sequenced isolates. This resulted in the removal of four wells leaving 68 in the analysis. For each of the 68 wells we then constructed parsimony trees manually to visualize the mutational relationships among the isolates.

Fitness assay

Because 15 of the 68 accepted wells contained one or more isolates where the ancestral mutation had reverted to wild type, we were led to question our original assumption that all background mutations were beneficial. We therefore conducted a follow-up set of competition fitness assays under conditions that matched the adaptation experiment. In these, we competed each background mutation (plus another mutation that commonly arose, A–G at site 1910, which we denote 1910aG) against wild type ID11 for six passages and estimated fitness from changes in mutation frequency. For each mutation we seeded three wells with a 90:10 mutant to wild type ratio and, in three more wells, the reciprocal ratio of 10:90. Sampling was done at passages 0, 2, 4 and 6 and population sequencing was conducted using the techniques described above except primer concentrations were manipulated to render 24-fold greater coverage in the region of the ancestral mutation.

Statistical analyses

Time to lysis assay and analysis

The result from Method 3 in conjunction with evidence in the literature led us to hypothesize that a number of the mutations that arose in the experiment delayed lysis. To test this, we identified four putative lysis-delaying mutations (1910aG, 1911cT, 2131cT and 2134tC) where we had one or more backgrounds both with and without the individual mutation in our set of sequenced isolates. We then performed an assay for the time to lysis where phage were added to cells in shaking flasks at 37 °C, sampled at one minute intervals and plated to determine when titer begins to rise. We then averaged the estimated burst size over replicates for each time point for each genotype and calculated standard errors.

Quality control

Contamination was a concern in our experiment because adaptations involved 100 transfer events for each of three 48-well plates using a multichannel pipette. We flagged potential contamination events based on four criteria: (1) the expected background mutation was missing from one or more isolates from the well, (2) in those isolates, one of the other background mutations appeared, (3) the putative contamination isolates carried other mutations (besides the other background mutation) that could link them to another well on the same plate and (4) reversions were otherwise rare in wells with the background we expected to find in the well. Upon obtaining our sequence data, one well showed all four red flags and was removed from the analysis. Three other wells met criteria 1 and 2, but not 3 and 4. It was impossible to know if the results from these wells was real evolution or represented contamination. Given the overall scarcity of evidence for contamination we suspect these represent real evolution. We therefore present our analysis with the wells included, but we also reran the analysis without them to confirm that their inclusion or exclusion has no qualitative effect on our results and conclusions.

Overview

The goal of this research was to assess parallelism and epistasis during adaptive evolution. The study, however, did not unfold as expected and we therefore begin with an honest summary of what transpired. The experimental design was to begin with nine different first-step beneficial mutations, allow each to further adapt, and have enough replication to formally assess how the accrued mutations depend on initial background, de novo mutations and chance. The nine background mutations, obtained in previous work (Rokyta et al., 2005), all arose and were highly beneficial in flasks. To achieve the desired replication (8 replicates per background) we shifted from passaging in flasks to sealed microtiter plates and transferred a fixed volume instead of a target number of phage. Titering during the experiment revealed that initial low MOI conditions in the plates rapidly gave way to high MOI ones. Sequencing endpoint isolates then revealed that 15 of the 68 wells (4 wells were removed due to sequencing failures) contained reversions, leading us to question our assumption that the background mutations were beneficial. Subsequent fitness assays and analysis (see Supplemental Information) led us to conclude that these background mutations were adaptive in a low MOI environment but not in a high MOI environment.

We analyzed the data from the perspective of parallel evolution and two major patterns emerged. First, we found that wells converged on the same phenotypic change. Most isolates in most wells, regardless of background, had exactly one mutation in the D-promoter. We further found a small number of other mutations just downstream from the D-promoter (in the coding region of genes D and E) that were mutually exclusive of each other and of the D-promoter mutations. The two most frequently observed D-promoter mutations have also been observed in experimental evolution of the bacteriophage ϕX174 (Wichman et al., 1999; Wichman, Millstein & Bull, 2005) and shown to down-regulate genes D and E (Brown et al., 2010). This led us to hypothesize that our D-promoter mutations and the downstream changes in the D/E genes were down-regulating expression or disrupting function of the lysis protein E and thereby delaying lysis. We reasoned that delay of lysis might be favored in the high MOI environment in which we passaged because it could increase burst size with little cost (since subsequent infections are unlikely at high MOI). We conducted an assay with several of the putative E-down-regulating mutations on several different genetic backgrounds and confirmed that they do indeed delay lysis. To quote Stephen Stills (who quoted Billy Preston), if you can’t be with the one you love, love the one you’re with (Stills, 1970).

Second, we found that the prevalence of parallelism depends on how and at what biological level we measure it at. From the population perspective, the proportion of nucleotide mutations shared between independent wells was 28%. If we redefine shared changes to be mutations in the same codon or regulatory feature, this value grows to 38%. Alternatively, we can look at individual mutations and ask how often they occur in more than one well. Overall, 36 of the 112 mutations (32%) were observed in more than one well. This moderate value, however, obscures how uneven the distribution of occurence was: 76 of the 112 mutations (68%) were found uniquely in one well and another 27 (24%) were found in between two and four wells. By contrast, 62 of the 68 wells (91%) we included in the analysis had a mutation in the D-promoter. We show below that the D-promoter mutations as well as a number of mutations just downstream of it delay lysis. Thus, at the phenotypic level, at least 64 of the 68 wells (≥94%) had mutations that delay lysis. Depending on the level of biological organization considered, parallelism might be characterized as low as 28% or over 90%, highlighting the need be attentive to the many biological levels at which the process occurs.

Mutational dynamics

After quality control (see Supplemental Information) we had 338 sequenced isolates distributed among 68 wells (4.97 isolates per well). Wells on average contained 4.5 mutations (range 2–9, sd = 1.5). We constructed parsimony trees for each well, in part to visualize the data, and in part to assess the types of mutational dynamics within the wells. The trees indicate that complex interference dynamics prevailed during the experiment—as we would expect in large populations—where multiple mutations of similar effect arise, compete and rarely fix before more mutations arise and begin increasing in frequency. Figure 1 illustrates the basic phylogenetic patterns and their prevalence. In just 4 of the 68 wells (6%), all of the sampled isolates were identical; this is consistent with, though not exclusive to, a simple model of sequential selective sweeps (Gillespie, 1984) (Fig. 1A). In the remaining 64 wells (94%) the population was not fixed. Models of adaptation differ with respect to the background(s) upon which secondary beneficial mutations arise. In the classic clonal interference model of Gerrish & Lenski (1998), multiple single mutations may arise and compete on a previously fixed background, but one of these secondary mutations fixes before the process repeats. In our data, 18 of the wells (26%) are consistent with this model (i.e., observations are confined to single tip mutations and their shared background; Fig. 1B). The other wells showed more complex, multiple-mutation dynamics (Desai & Fisher, 2007) where beneficial mutations arise on unfixed backgrounds. In 5 of the wells (7%) the trees were unbranched, consistent with mutations arising sequentially on the single most-fit background (Fig. 1C) while in 41 of the wells (60%) the tree was branched in ways that suggest competing multi-step lineages (Fig. 1D) (Miller, Joyce & Wichman, 2011; Desai, Fisher & Murray, 2007). The entire set of phylogenies are presented in the Supplemental Information. These findings add to a growing body of evidence (e.g., Kao & Sherlock, 2008; Miller, Joyce & Wichman, 2011; Desai, Fisher & Murray, 2007; Lee & Marx, 2013; Lang et al., 2013; Frenkel, Good & Desai, 2014) that mutational dynamics in large populations usually involve complex interference dynamics.

Characterizing the mutations

A total of 112 de novo mutations were observed across all wells. These spanned most of the genome and fell within all but one gene. However, mutations were highly uneven both in where and how many independent times (wells) they occurred (Fig. 2). By both of these measures, activity was concentrated in the D-promoter and in genes D and E (which have overlapping reading frames). It was also moderately high in the J terminator and in gene F. Of the 112 mutations, 15 were in regulatory regions or known to have regulatory effects, 76 were non-regulatory and non-synonymous, and the remaining 21 were silent. Multiplying each mutation by the number of wells it occurred in, we have 281 de novo mutational events in the experiment; 95 (33%) of these were mutational events in the D-promoter, 80 (28%) were in genes D and E, 57 (20%) were in gene F and 49 (19%) of the events were elsewhere in the genome. Matrix files with all the information about every observed mutations in every isolate from every well and the ancestral state of that well are available as Supplemental Information.

Neutral variation will tend to reduce parallel evolution. This leads us to ask, how much of the observed variation is adaptive vs. how much is neutral? To answer this, we note that in our data there are two signatures that a mutation is beneficial: (1) when it appears in multiple wells independently, and (2) when it rises to moderate or high frequency unaccompanied by another mutation. The phylogenies from each well contain information about which mutations meet this second criteria as explained in the Supplemental Information (subsection Identifying Neutral vs. Beneficial Mutations). Based on these lines of evidence, at least 46 of the 112 mutations (41%) are adaptive (Fig. 3A). But because many of these changes occur in multiple wells, the 46 adaptive mutations comprise 215 of the 281 de novo mutational events observed in the experiment. Thus, at least 77% of the mutational events are adaptive (Fig. 3B). The Venn diagram also partitions the mutation and event counts into silent changes on the one side vs. those that are non-synonymous and in regulatory regions (e.g., in the D-promoter) on the other. Two things jump out. First, 41% of the mutations and 17% of the mutational events are both nonsynonymous and lack any evidence of being adaptive. Second, there are only 22 silent mutations in the data. Three of these, corresponding to 16 (or 6%) of mutational events, are adaptive.

Parallel evolution

The focus of this study is to assess the level of parallel evolution. We begin there, digress into one of its underlying determinants–epistasis–and then follow this topic back around to parallel evolution. Parallel evolution refers to the same change occurring in independent trials of evolution. Although conceptually simple, quantifying it involves some nuance as illustrated by the following three questions. First, what constitutes the same change? We might define shared changes at several biological levels: matching nucleotide change, matching amino acid changes, a change at the same residue, in the same protein domain, in the same gene, in the same regulatory element, part of the same regulatory network, or affecting the same phenotype. Second, should we quantify the frequency of parallelism from the perspective of the mutations or of the replicate populations? In other words, we can ask how often change x co-occurs in two trials, or we can ask of a random pair of replicates, how many or what proportion of the changes are shared. Third, how similar must the starting points of the trials be to consider them replicates? One thing we wish to emphasize here is how different the answers can be.

Parallel evolution by mutation

Of the 112 DNA mutations observed in the experiment, 36 (32%) were observed in more than one well whereas 76 (68%) were observed in only one well. Six of the observed mutations were reversions. We begin by focusing on the non-reversions as this is where the opportunity for evolution is the same across wells. When we plot the distribution of well occurrences by mutation, we find it is highly uneven. Two D-promoter mutations are found often (in 53 and 27 wells) while 103 mutations (92%) are found in four or fewer wells (Fig. 4A). With an observed well frequency of 0.78 (53/68), the probability of finding the D-promoter mutation 1910aG in each of two randomly selected wells is 0.61 (0.78²). By contrast, a mutation like 2397aG that is found 8 times (frequency 0.12) has just a 0.014 probability of co-occurring in a randomly selected pair of wells. When we do this calculation for all 106 mutations and take an unweighted average, we find that a randomly selected mutation in this experiment has only a 0.009 probability of occurring in each of two randomly selected wells.

Although the D-promoter mutations 1910aG and 1911cT arise far more often than other mutations, the distribution in Fig. 4A probably understates the level of parallelism. If we assume that all changes in the D-promoter have similar regulatory effects (an assumption we return to later), then it is appropriate to group them together as a single mutational feature. Upon doing so we find that 62 of the 68 wells (91%) have a D-promoter mutation and that the probability of two random wells each having one is 0.83. We also clustered all DNA mutations that alter the same codon, but because very few mutations affect the same residue, this had virtually no impact on the distribution of counts.

We now turn to the reversions. Each background had between 6 and 8 replicate opportunities to revert in the experiment. The number of observed reversions differed dramatically by background, ranging from zero in three backgrounds up to seven in background F416 (Fig. 4B). A likelihood ratio test confirms that background has a significant effect on the probability of a reversion (p = 0.0002). Thus, for some backgrounds, reversions were a major contributor to parallel evolution and at others they were not. This finding also highlights that genetic background may influence how likely two wells are to share a mutation. We provide to a more detailed analysis of these mutational interactions in the epistasis section below. More immediately, we use this result to justify why, in the next analysis on parallelism between wells, much of our focus is within (not between) backgrounds.

Parallel evolution by well

In the last section we asked how likely a given mutation was to appear in two independent wells. We now change the conditioning and ask, for two independent wells evolved from the same background, what proportion of the mutations observed in them are shared? The short answer, as we will show in the next paragraph, is 28% at the nucleotide level and 38% when mutations are defined at the regulatory/codon level. At the DNA level, this is somewhat lower than the 0.5 value observed in a ϕX174 experiment (Wichman et al., 1999), but greater than 2–3% values observed in E. coli experiments (Woods et al., 2006; Tenaillon et al., 2012). The increased level of parallelism at the regulatory/codon level is qualitatively similar, if less dramatic, than the findings in E. coli where parallelism increases an order of magnitude at the gene/operon level. (Because there are only 11 genes in the bacteriophage, it is not all that meaningful to ask how often lines share changes in the same gene.)

To make our analysis more precise, we define between-well parallelism as the probability that a mutation observed in one well was also be found in a second well. We calculated observed parallelism for all pairs of wells initiated from the same background, denoted i and j, by determining the proportion of mutations from j found in i, the proportion of mutations from i found in j, and averaging. At the nucleotide level, mean within-background parallelism ranged from 0.15 in F421 to 0.46 in F416, with an average of all nine backgrounds of 0.28 (SE = 0.03; Fig. 5A). The reason F416 shows such high parallelism is that reversions occurred in 7 of 8 cases (Fig. 4B) and, in each case, mutation 2534gT also arose. No other background experienced more than 3 reversions. When reversions are removed from the dataset, parallelism at F416 drops to 0.38 and the global mean drops just slightly to 0.27. We conducted a randomization test that accounts for the nonindependence of pairwise comparisons and found that, with reversions excluded, there is no evidence that parallelism differs among backgrounds (p ≈ 0.58).

We have already shown that the D-promoter mutations are found in most wells (Fig. 4A). To quantify their contribution to parallelism, we removed the D-promoter sites and reanalyzed the data. Within background parallelism falls from an average of 0.28 to 0.13. Thus, slightly more than half of the nucleotide level parallelism comes from the D-promoter mutations.

A broader view of parallelism is to categorize mutations as the same when they occur in the same regulatory element or codon. When we do this, parallelism increases from an overall mean of 0.28 to 0.38 (Fig. 5B). The randomization test indicates that F416 is different from the other backgrounds (p ≈ 0.013), but the other eight backgrounds show no evidence of having different levels of regulatory-level parallelism (p ≈ 0.228). Removing reversions from the dataset has a minimal effect on mean parallelism (0.378 falls to 0.376). When F416 is removed (but reversions are left in), mean parallelism falls from 0.38 to 0.35. When the D-promoter mutations are again excluded from the analysis, the overall mean of this broad parallelism falls to 0.17. We regard 38% to be our single best characterization of parallelism in this experiment, with slightly more than half of this owing to parallel changes in the D-promoter.

Epistasis

Epistasis refers to the way in which the genetic background influences a mutation’s fitness effect. Parallelism is increased when beneficial mutations have approximately the same effect on all genetic backgrounds and it will decrease when the fitness effects (including their sign) depend strongly on the genetic background. By initiating our replicate adaptations from nine different backgrounds, our experiment was explicitly designed to test how epistasis influences where mutations arise and, thereby, the amount of parallel evolution in our system.

Effect of Background

We have already seen that reversion rate is significantly different among backgrounds (Fig. 4B). In this section, we begin looking for background epistasis by asking whether there is a difference in parallelism within- vs. between-backgrounds once reversion are removed. Figure 6A shows a heatmap of all pairwise well comparisons at both nucleotide (red) and regulatory/codon (blue) levels. Reversions have been removed in this analysis. Notice that visually the within-background comparisons along the diagonal are not much different than those off the diaganol by either metric of parallelism. This suggests that parallelism is similar within and between backgrounds. However, when we pool all observations within vs. between, we find that the distribution of within-background parallelism is slightly greater than between (Fig. 6B). Thus, parallel mutations are slightly more likely to occur within wells of the same background than of different backgrounds. A formal randomization test accounting for nonindependence of the observations confirms that the difference in means is significant (p < 0.01). The implication of the within- vs. between-background analysis of Fig. 6 is that while the occurrence of mutations is highly stochastic, genetic background does exert a slight but significant influence over where mutations appear.

A second, more powerful, way to assess the role of background epistasis is to ask whether background affects the individual probability of a mutation arising. More formally, does a model where the probability of a mutation occurring depends on background fit the data significantly better than a null model where background does not matter? Given the parallelism results just presented, it is not surprising that in a likelihood ratio test (LRT) of these competing hypotheses the null can be rejected (p < 0.001). The main advantage of the LRT approach is that the disparity between the alternative and null hypotheses (i.e., the difference in the log-likelihoods, or ΔlnL), can be broken down by mutation. Figure 7 shows the number of times each mutation appears on each background and, to the right, ΔlnL for each mutation. The larger bars are the mutations most responsible for the association of mutation with background. We removed mutations in a cumulative manner starting with the largest ΔlnL value and reran the LRT until we failed to reject. This revealed that just four of the mutations (red bars in the figure) are principally responsible for the significant relationship: 2534gT, 2397aG, 2332aG, and 2100tC. Thus, we again find that background epistasis exists, but the effects are substantial for only a small minority of the mutations.

Epistasis between de novo mutations

Another type of epistasis involves combinations of new mutations. For our data, the signal for this is when mutations appear together significantly more or less often then we would expect by chance. We conducted an analysis on mutation pairs, asking if they tend to co-occur or avoid each other. The results, summarized in Fig. 8, revealed two interesting patterns. First, despite finding D-promoter mutations in 62 of 68 (91%) wells and in 273 of 338 (81%) isolates, no two D-promoter mutation are ever found together in the same isolate (Fig. 8). The significance of this is seen clearly in the figure by focusing on the D-pro block in the bottom row which represents the combined set of all the D-promoter mutations. A bootstrap randomization indicated we expect D-promoter mutations to co-occur 19 times; observing 0 is highly nonrandom (p < 0.01). Second, several other mutations are never found with a D-promoter mutation. Most notably, mutation 2134tC is found in 11 different wells and is expected to co-occur with a D-promoter mutation 7 times, but it never does. Mutations 2094aG, 2131cT and 2158aG are also never observed with a D-promoter mutation, but these mutations are rare enough that the observed repulsion is not significant. Finally, the analysis shows several significantly positive associations that suggest synergistic epistasis: 2534gT co-occurs with D-promoter mutations and with 2630gT, and 2577cT occurs with the D-promoter mutation 1910aG more often than expected by chance.

Parallel evolution at the phenotypic level

The most striking result from analysis of epistasis is that while D-promoter mutations are pervasive, they never co-occur in the same isolate, and a number of other mutations are also mutually exclusive of the D-promoter mutations and each other. One parsimonious explanation of this pattern is that these mutually exclusive mutations have closely related phenotypic effects and that having two such mutations is either neutral or deleterious. We now show that the phenotype involves delay of lysis. We postulate that selection favored lysis delay because, under high MOI conditions, lysing early has a diminished benefit (finding another host cell is improbable) compared to the advantage of exploiting the current host more effectively, i.e., producing more progeny.

The first evidence of what selection is doing here comes from experimental evolution on ϕX174. In chemostat adaptation experiments where high MOI conditions prevail, D-promoter mutations are commonly observed (Wichman et al., 1999; Wichman, Millstein & Bull, 2005). Brown et al. (2010) studied four D-promoter mutations and found that they all reduce gene D/E mRNA transcript levels on the order of 75–80%. Two of the mutations studied there are the same as the two most commonly observed in our data: 1910aG and 1911cT. Since the lysis protein E is encoded in this transcript, protein E levels are presumably reduced by these mutations. We therefore suspected that our D-promoter mutations are down-regulating the lysis protein in our experiment and thereby delaying lysis.

This led us to examine the 65 isolates that lack a D-promoter mutations. Might these be achieving similar ends by different means? We found that 56 of these 65 isolates (86%) have mutations that could reasonably affect the E protein in that the mutations (i) are either just upstream of the E protein coding region or in the gene itself and (ii) are never found in an isolate with a D-promoter mutation (46 of 56 have the aforementioned mutations 2094aG (5 isolates), 2131cT (4 isolates), 2134cT (26 isolates) and 2158aG (11 isolates) while the other 10 have mutations found in just one well: 2147tA (2 isolates), 2335gT (1 isolate), 2361cG (5 isolates), 2372gT (1 isolate) and 2407aT (1 isolate)). Thus, only 9 of the total 338 isolates (2.7%) had neither a mutation in the D-promoter nor a change that might reasonably effect the E protein.

These observations led us to do lysis assays on a subset of mutations we hypothesized were delaying lysis: 1910aG and 1911cT (the two common D-promoter mutations found in 182 and 62 isolates respectively), 2131cT (found in 4 isolates) and 2134cT (found in 26 isolates). These mutations were selected in part because of their abundance (representing 274 of 338 isolates, or 81%) and because we already had sequenced isolates with and without each mutation. For three of them (all but 2131cT) we had isolates allowing the with/without comparison to be made on multiple backgrounds. This let us assess phenotypic effect of each mutation across multiple genetic backgrounds. The results, presented in Fig. 9, confirm our hypothesis. The D-promoter mutations delay lysis 5–8 min depending on the mutation and the background. Mutations 2131cT and 2134cT delay lysis on the order of 12–13 min. The patterns are very similar across background suggesting that this phenotypic effect is not very background dependent. Thus, we have strong evidence that over 80% of the isolates are delaying lysis and, based on the arguments above, there is reason to suspect that most of the remaining isolates have delayed lysis as well. At the well level, at least 64 of the 68 wells (≥94%) have mutations that delay lysis or are strongly suspected of doing so.

Given this dramatic delay in lysis time, we were interested in what molecular mechanism(s) underlie the delay. With the D-promoter mutations, we hypothesized that the mutations change the binding affinity of the RNA polymerase which alters expression of the lysis protein, protein E (Brewster, Jones & Phillips, 2012; Kinney et al., 2010). For mutations downstream of the D-promoter within the D/E genes, we hypothesized that delay is driven by other means of delaying protein E expression. We considered the possibility that E is down-regulated by mutation to rare codons (Sharp, Emery & Zeng, 2010), ribosomal pausing at mutationally introduced Shine-Delgarno-like sequences (Li, Oh & Weissman, 2012) or mutationally altered transcript stability (Agashe et al., 2013). As we detail in the Supplemental Information, our analyses failed to find strong support for any of these hypotheses by themselves. We continue to suspect that lysis delay is driven by change in protein E expression. We speculate that our inability to uncover drivers may come from having relatively few mutations, from using models (as opposed to molecular data) for exploring the mechanisms and the possibility that multiple mechanisms may be involved across and even within individual mutations.