A new practical guide to the Luria–Delbrück protocol

Mutation Research 781 (2015) 7–13

Contents lists available at ScienceDirect

Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis journal homepage: www.elsevier.com/locate/molmut Community address: www.elsevier.com/locate/mutres

A new practical guide to the Luria–Delbrück protocol Qi Zheng Department of Epidemiology and Biostatistics, Texas A&M School of Public Health, College Station, TX 77843, United States

a r t i c l e

i n f o

Article history: Received 17 June 2015 Received in revised form 25 July 2015 Accepted 23 August 2015 Available online 28 August 2015 Keywords: Plating efficiency Relative fitness Likelihood ratio test Mutation rate

a b s t r a c t Since 2000 several review papers have been published about the analysis of experimental data obtained using the Luria–Delbrück protocol. These timely papers cleared much of the confusion surrounding various methods for estimating or comparing mutation rates. As a result, today the fluctuation test is more widely applied with much improved accuracy. The present paper provides guidelines on a few remaining problems that continue to baffle mutation researchers. Among the issues addressed are incomplete plating, relative fitness, and comparison of experiments where average final cell population sizes differ. It also offers a fresh view on the estimation methods that are based on the sample median. © 2015 Elsevier B.V. All rights reserved.

1. Introduction One Saturday evening in early 1943, having mulled for several months over the then ongoing debate about whether mutations were spontaneous or adaptive, Salvador E. Luria suddenly found an experimental method to solve the fundamental controversy. After Luria completed the experiments, Max Delbrück devised a mathematical framework within which Luria’s experimental data were analyzed and interpreted. Thus was born the well-known Luria–Delblrück (henceforth LD) protocol. The experiments of Luria came to be known as the fluctuation tests or fluctuation experiments [26]. From the classic paper that reported the protocol [25], it is immediately clearly that the LD protocol not only sheds new light on the controversy about the origin of mutation, but also offers a novel approach to the determination of microbial mutation rates under laboratory conditions. While the exact role that the LD protocol played in the controversy is still debatable today, a consensus had been reached decades ago that the LD protocol is the method of choice for estimating and comparing microbial mutation rates in the laboratory. Under the LD protocol, the experimentalist seeds a number of tubes of liquid culture with small numbers of wild-type cells. After incubating the liquid cultures, the experimentalist transfers (plates) either the entirety or a part of the contents of each tube onto a dish of solid, selective culture. Because the wile-type cells are rapidly eliminated by the selective culture, the colonies that appear during the second incubation period represent the mutants

E-mail address: [email protected] http://dx.doi.org/10.1016/j.mrfmmm.2015.08.005 0027-5107/© 2015 Elsevier B.V. All rights reserved.

that exist immediately prior to plating. In essence, the LD protocol renders the invisible mutations partially visible to the experimentalist, which in turn makes possible the estimation and comparison of microbial mutation rates via mathematical reasoning. Using the LD protocol to study mutation rates seems more straightforward than using it to solve a fundamental controversy. However, even for this less ambitious purpose, the LD protocol has presented its users with peculiar conceptual and computational difficulties since its inception. At times, even the most diligent users of the LD protocol were baffled by these difficulties. For example, Newcombe [30] designed sophisticated experiments to help understand the discrepancies between estimates obtained by the p0 method and those by the method of the mean, the two major methods available at the time for determining microbial mutation rates using the LD protocol. But today the quest to fathom the method of the mean still seems unfinished. It was amid such persistent confusion that the first guide to the LD protocol appeared in early 1988 [20]. Mutation researchers at that time still relied exclusively on the p0 method and the method of the mean, which had been in use for some 45 years since Luria and Delbrück applied them in the world’s earliest fluctuation experiments. (There were sporadic methodology developments, e.g., Ref. [21], but they were not accessible to most bench scientists.) The thought-provoking prose of Kendal and Frost [20] was not so much a practical guide to what one could do as a cautionary note against what one should not do. The means by which the experimentalist can avoid the pitfalls identified by Kendal and Frost did not exist yet. However, in the same year a far-reaching event in mutation research occurred that would eventually create a genuine need

8

Q. Zheng / Mutation Research 781 (2015) 7–13

for a practical guide to the LD protocol. The publication of Ref. [6] not only revived the spontaneous-adaptive mutation debate that was directly responsible for the coming into being of the LD protocol, but also spurred a surge of research into new methods (a method of point estimation is commonly known as an estimator) for estimating and comparing mutation rates using the LD protocol. Among the important contributions to methodology development during this “renaissance” period are Refs. [1,19,27,38,39]. The appearance of many new methods posed a great, albeit welcome, challenge to practitioners. The timely review by Rosche and Foster [35] played an important role in the assessment of the estimation methods available at the end of the 20th century. In particular, this review greatly increased awareness of the optimal properties of the maximum likelihood (ML) method, thereby gradually weaning bench scientists from a then popular yet unreliable method—Delbrück’s method of the mean. A second review by Foster [17] complemented the first one by suggesting additional methods for computing confidence intervals (CIs). Recognizing the importance of the LD protocol in the study of drug-resistant microbes, Pope et al. [33] reinforced messages from the two reviews by condensing them into a minireview. Surprisingly, the distillation process discarded the ML method, putting undue emphasis on methods based on the sample median. In 2009, a new distillation process occurred. The web-based tool FALCOR [15] by Hall et al. [18] made the ML method (along with the Lea–Coulson (LC) method) easily accessible via the Internet, occasioning a second surge of interest in the LD protocol in genetic research. A spate of applications of the LD protocol were facilitated by this convenient tool, e.g., Refs. [2,7,9,16,22,24,29,34,37]. Finally, a large-scale simulation study by Couce and Balzquez [8] corroborated the optimal properties of the ML method, further increasing the popularity of the ML method and discouraging the use of two median-based methods. These five papers, although spanning more than a decade, acted collectively as a guide for practitioners in the use of estimation methods available at the end of the 20th century. What, then, is the point of another practical guide? The motivation for the present guide was the observation that a considerable portion of the recent applications of the LD protocol still suffers from the use of incorrect or suboptimal methods. This situation was caused by two facts. First, some existing methods, e.g., Drake’s formula [12], were not fully understood at the end of the 20th century. Second, modern methods for solving some routine problems did not exist at the end of the 20th century. Among such “routine” problems are partial plating, fitness cost and comparison of mutation rates. Most of these problems can be solved relatively easily when the likelihood paradigm is fully adopted. Although use of the likelihood approach in the context of the LD protocol can be traced to as early as Lea and Coulson [23], systematic use of the likelihood approach in computing point and interval estimates and performing statistical tests was a relatively recent development [45,47,50]. The mutation rate can be directly tackled using, e.g., the Haldane model [46]. However, it is historically more logical and mathematically less cumbersome to focus on m, the expected number of mutations per culture. Thus, this paper focuses exclusively on the LC model [23] and a simple variant. The mutation rate is simply m divided by Nt , the final cell population size. Due to historical reasons, an estimate of m was often converted to a mutation rate by log(2) × m/Nt . The log(2) factor was originally introduced by mistake [44] and is no longer necessary (in mathematics, the natural logarithm is commonly denoted by the symbol log, in stead of ln; the present paper follows this convention). Before recommending new methods, I explain why an old approach should be regarded as obsolete, reinforcing an observation by Couce and Balzquez [8].

2. Demystifying the concept of the likely average In a recent investigation, Oide et al. [31] chose two medianbased methods to estimate mutation rates, according to whether or not the parameter m was believed to be greater than 30. These two methods are commonly known as the Lea–Coulson method [23] and Drake’s method [12], respectively. A lack of understanding of the close relation of these two methods to Delbrück’s method of the mean was responsible for the use of this unsound rule (in view of a footnote in the original paper [25] and a later recount by Luria [26, p. 79], it seems safe to assume that the idea of the method of the mean was Delbrück’s alone). Because the method of the mean relies on Delbrück’s concept of the likely average, I first give a description of this concept in present-day terminology. First, focus on one test tube. Let N(t) = N0 eˇt be the number of wild-type cells at time t, and let X(t) be the number of mutants at time t. Viewing the number of mutants X(t) as a shot noise process [42], Delbrück concluded (e.g., via Campbell’s theorem) that E[X(t)] = uN0 teˇt , where u is the probability of mutation per cell per unit time. It thus follows that E[X(t2 )] E[X(t1 )] − = u(t2 − t1 ). N(t2 ) N(t1 )

(1)

Equating the two expected mutant frequencies in the above equation with the corresponding observed mutant frequencies f1 and f2 , and noting that t2 − t1 = ˇ−1 log(N(t2 )/N(t1 )), one has p≡

u f2 − f1 = . ˇ log(N(t2 )/N(t1 ))

(2)

Here p is the probability of mutation per cell division, or the mutation rate [42]. Because X(t) in a single test tube cannot be measured twice, Delbrück considered measuring the first mutant frequency f1 “in theory” at an intuitively appealing time t1 . Now, allow the experiment to have C test tubes. Delbrück chose t1 = −ˇ−1 log(pCN0 ), to achieve pCN(t1 ) = 1. Delbrück further postulated that at this critical point, E[X(t1 )] = 0. Consequently, f1 = 0. Thus, with this fixed t1 and an arbitrary t = t2 > t1 , Eq. (2) reduces to p=

f2 . log(pCNt )

(3)

If E[X(t)] = r, then f2 = r/Nt . Therefore, Eq. (3) becomes the celebrated Eq. (8) of the classic paper [25]: r = pNt log(Nt Cp).

(4)

For numerical convenience, Eq. (4) is often recast in term of m = pNt : m=

r . log(Cm)

(5)

Because the sample mean of a mutant distribution can be unreliable, the sample median was soon widely accepted as a substitute for the sample mean r in Eq. (5). This seemingly effective remedy diverted attention from another flaw of Delbrück’s method—the concept of the likely average. Newcombe [30, p. 453], calling the special time t1 “an arbitrary time,” was perhaps the earliest to question the usefulness of the above reasoning. However, it is more instructive to examine the arbitrary choice of zero as the expected number of mutants at time t1 . One can argue that setting E[X(t1 )] = 1 would be equally appealing, because at this epoch a mutation is likely to have occurred but the first-generation mutant has not yet generated its offspring. Thinking along this line, one might consider E[X(t1 )] = b for an arbitrary, small number b. This reasoning, coupled with the identity

Q. Zheng / Mutation Research 781 (2015) 7–13

9

b=1.75

b=2.0

ML

the plating efficiency. Partial plating was practiced in 9 out of the world’s 10 earliest fluctuation experiments [25]. Today, partial plating remains a popular option. Plating efficiencies of 0.1 and 0.2 are the most favored ones. Prior to the publication of Ref. [47], two empirical methods were commonly employed to account for partial plating. The first is the method due to Jones et al. [19]: ˆ 0.5 / − log 2 , ˆ log(0.5 /) − log(log 2)

b=1.5

ˆ = m

b=0.0

b=0.5

b=1.24

where ˆ 0.5 denotes the median of the mutant count data (note that this method is not applicable when ˆ 0.5 = 0). The second method, suggested by Stewart et al. [38], was articulated in Ref. [35], and further popularized by FALCOR’s website [15]. This method requires two steps. First a temporary estimate of m is obtained by assuming that plating is complete ( = 100%). The resulting (ML) estimate is denoted by mobs . Second, an estimate of a so-called “actual” m is derived by applying Eq. (41) of Stewart et al. [38]: 20

30

40

50

60

70

mact =

80

Fig. 1. Ten thousand 15-tube fluctuation experiments were simulated with m = 30.0. Six Ab estimators and the ML estimator were applied to the simulated data. The distributions of the resulting estimates of m are summarized by boxplots. The A1.24 estimator is commonly known as the Lea–Coulson estimator, and the A0 estimator is also called Drake’s estimator. The vertical dashed line indicates the true parameter value m = 30.0.

−1 m .  log() obs

(9)

However, it has so far escaped attention that Eq. (9) does not hold exactly. This is most easily understood by examining the two relevant probability generating functions (PGFs). The PGF of the distribution of the actual number of mutants is of the form [47]



G1 (z; m, ) = exp m pCN(t1 ) = 1, leads to the “imaginary” mutant frequency at time t1 being expressible as f1 = b/N(t1 ) = pbC. Eq. (3) in turn becomes p=

f2 − pbC . log(pCN t )

(6)

In the same vein, Eq. (6) can be rewritten as m=

(8)

r − bCm , log(Cm)

(7)

which generalizes Eq. (5). Here r is an estimate of E[(X(t)], which can be either the mean or the median number of mutants. The LC method [23], m log(m) + 1.24m = r, along with the Drake formula m log(m) = r [12], suggests that there might be an optimal value of b for each given C. That is, C can be regarded as a parameter independent of the actual number of cultures. For convenience, let Ab denote the estimator (7) with C = 1 and r being the sample median. Thus, A1.24 is the LC estimator, and A0 is the Drake estimator. A simulation (summarized in Fig. 1) indicates that 1.24 is close to the optimal value of b for C = 1. It is also clear from the simulation that the A0 estimator is noticeably biased and has a larger variance than the ML estimator. These problems of the A0 estimator do not seem to diminish when m exceeds 30. For an arbitrary integer C > 1, the behavior of the estimator (7) is unclear. With the ML method readily available, it is no longer wise to use these methods that are based on the sample median [8].

(1 − z) log[(1 − z)] 1 + z



,

(10)

where =

 1−

.

The validity of Eq. (9) rests on the assumption that the data can also be regarded as generated by an LC distribution. This imagined distribution has the PGF of the form G2 (z; m∗ ) = exp



m∗

1 z





− 1 log(1 − z)

.

(11)

If the claim were true, one could chose a parameter value of m in G1 such that the two generating functions coincide. This is not possible, considering just the first two terms in the series expansions of the above two generating functions. First, for the probability of zero mutants, the two probabilities are respectively P0∗ = exp(−m∗ ) and

P0 = exp

 m log   1−

(12)

.

To match the above two terms, one has to chose m = ( − 1)m*/( log ), which is Eq. (9). But the next two probabilities are P1∗ = P1 =

m∗ exp(−m∗ ), 2 m( − 1 − log ) (1 − )2

 exp

1+

m 1−





(13)

log  .

Substituting Eq. (9) into the second equation, one has 1 −  + log  , (1 − ) log 

3. Accounting for partial plating

P1 = m∗ exp(−m∗ )

Unlike somatic cells, bacterial cells are less fragile and their death during the incubation period and the ensuing plating process is generally regarded as negligible [20]. Therefore, imperfect (<100%) plating efficiency is most often caused by intentional partial plating. As Boe et al. [5] reported, plates with more than 512 mutants cannot be counted with precision (see also Ref. [1, p. 314]). The experimentalist is often compelled to plate only a portion (denoted  < 1) of the contents of each test tube. Here  is called

which does not agree with P1∗ in the first equation. When  is close to 1.0, Eq. (9) can be a useful approximation. For values of  that are common in practice, e.g., 0.1, the approach based on Eq. (9) can lead to noticeably biased estimates of m. A more confusing problem is the construction of CIs. The literature has been vague about this important issue. Most authors seem to rely on an empirical method that was originally proposed for cases where plating efficiency is perfect ( = 1.0). This method, based on

(14)

10

Q. Zheng / Mutation Research 781 (2015) 7–13

Table 1 Two methods of accounting for plating efficiency.



0.1

0.2

0.3

0.4

0.5

0.6

0.8

Median Coverage Length

6.6/8.0 73.4/94.7 4.5/7.0

6.7/8.1 78.0/94.9 4.5/6.1

6.9/8.1 83.0/95.0 4.6/5.7

7.1/8.1 87.1/94.8 4.7/5.5

7.3/8.1 90.6/95.3 4.8/5.4

7.4/8.0 92.5/94.8 4.8/5.3

7.7/8.1 94.3/94.9 5.0/5.1

Comparison of the two available methods of adjusting for plating efficiency: the method based on Stewart’s empirical formula (9) and a likelihood-based method. For each of 7 given plating efficiencies, 104 15-tube fluctuation experiments were simulated with a common value of m = 8.0. Three comparison indices were calculated: “median” refers to the median of the point estimates of the parameter m; “coverage” refers to the percent coverage of the true parameter value m = 8.0 by 95% confidence intervals; “length” refers to the median length of the 95% confidence intervals. Each entry is in the form of a/b, where a is a result obtained using the Stewart empirical approach, and where b is the corresponding result obtained using the likelihood approach described in Ref. [47].

CL95+ = log(m) + 1.96e−1.35×1.96 CL95− = log(m) − 1.96e1.35×1.96

(15)

Jones

an observation of Stewart [39], was refined by Rosche and Foster [17,35]. It starts with setting up a 95% CI for log(m), the upper and lower confidence limits being given by

1.225m−0.315 . √ n

(16)

Here, n is the number of cultures (sample size). Most authors seemed to extend the above approach to the case of incomplete plating by replacing m in Eq. (15) with mact obtained from Eq. (9). CIs thus obtained were tacitly accepted in the literature, but their statistical properties are unknown. The empirical estimate mact based on Eq. (9) is often mistakenly referred to as an ML estimate in the literature. A bona fide ML estimate of m can be calculated using the Newton–Raphson technique described in Ref. [47]. CIs for m can also be constructed using a classic likelihood-based approach [47]. To assess the empirical method, I performed the following simulation. For each of 7 prescribed values of plating efficiency ( = 0.1(0.1)0.6, and 0.8), a set of 104 15-tube fluctuation experiments were simulated. Other parameters were the same for the seven sets of experiments: N0 = 10, Nt = 2.0 × 108 and a mutation rate of 4 × 10−8 . Therefore, the true value of m is 8.0. Point estimates of m and CIs are calculated using both the empirical method and the likelihood method. Results are summarized in Table 1. As Table 1 indicates, because the empirical approach ignores the variability due to sampling, the CIs given by the empirical approach are narrower than the corresponding CIs given by the likelihood-based approach. As a consequence, the 95% CIs given by the empirical approach do not cover the true parameter value (m = 8.0) with the intended probability of 0.95. For example, in the case  = 0.1, the coverage rate is only 73.4%. As far as point estimates are concerned, the empirical approach tends to underestimate the true value of m. Fig. 2 depicts the distributions of estimates produced by the three available methods for the case  = 0.2. Both the Jones method and the Stewart empirical method are noticeably biased. Note in passing that in the above and subsequent simulations, any simulated mutant counts greater than 5 × 104 are truncated to 5 × 104 . 4. Accounting for fitness It is generally believed that a mutation that confers resistance to an antibiotic renders the mutant cell carrying that mutation less competitive in an environment free of the antibiotic. An important manifestation of reduced competitiveness is diminished growth rates of the mutant. The concept of relative fitness measures the growth rate of mutants relative to that of nonmutants. In a fluctuation experiment, antibiotic resistant-mutants, generated spontaneously, do not come into contact with the antibiotic prior to plating. Because the LC model [23] assumes that mutants

Stewart

=

ML

with

4

6

8

10

12

14

16

18

Fig. 2. Ten thousand 15-tube fluctuation experiments were simulated with m = 8.0 and plating efficiency  = 0.20. The Stewart method (9), the Jones method (8) and the ML method were applied to the simulated data to estimate m (assuming  = 0.2 known). The distributions of the resulting estimates are summarized by boxplots. The vertical dashed line indicates the true parameter value m = 8.0.

and nonmutants grow at the same rate, the effects of reduced relative fitness on the estimation of mutation rates have been largely ignored in the past. Exceptions are the likelihood-based algorithms proposed in Ref. [45] to estimate jointly the mutation rate and relative fitness, and to construct CIs for both the mutation rate and relative fitness, using only fluctuation assay data. However, relative fitness is more reliably measured by a separate fitness assay. Fitness assays, also referred to as competition assays, are increasingly common in studies of microbial mutation rates [2,4,16,29,36,40]. Due to lack of suitable algorithms, these fitness assays, although often performed in tandem with fluctuation experiments, were not regarded as an integral apart in the estimation of mutation rates. The present section addresses this situation. To begin with, I articulate the way relative fitness is defined in the context of the LD protocol. There are several definitions of relative fitness [32]. The most appropriate and convenient one for fluctuation experiments is the ratio of two Malthusian parameters, which is widely used in experimental evolution [10]. In a fitness assay [13,36], two strains of bacteria of comparable densities are inoculated, say N01 and N02 . After a predetermined period of growth, the numbers of the two strains (Nf1 and Nf2 ) are determined. The relative fitness is calculated as



ˆ = log w

Nf2 N02





÷ log

Nf1 N01



.

(17)

Q. Zheng / Mutation Research 781 (2015) 7–13

11

ˇmutant . ˇwild-type

(18)

Here ˇmutant and ˇwild-type are respectively the growth rates of mutants and nonmutants. Definition (18) is equivalent to that given in Ref. [4]. Logarithmic (i.e., exponential) growth of both types of cells is tacitly assumed in this definition of relative fitness. Note that some authors [16] use the difference ˇmutant − ˇwild-type as relative fitness, which is not convenient in the context of a fluctuation experiment. To integrate information about w from a fitness assay into the estimation process, one first needs to modify the popular LC model [23]. This has already been accomplished by Mandelbrot [28], who developed the first stochastic model that extends the LC model to accommodate differential growth. Independently, Koch [21] devised largely empirical methods for inference under the same model. I thus call the model first studied by Mandelbrot the Mandelbrot–Koch (MK) model. Under the MK model, the PGF of the exact mutant distribution is (e.g., Eq. (23) in Ref. [43])

 G(z; m, w, ) = exp

m

−m + Bω (k, 1 + w−1 )z k w ∞

,

(19)

k=1 w

ω = 1 − (1 − ) ,  = 1 − N0 /Nt and Bz (a, b) = b−1 t (1 − t) dt is the incomplete beta function. For appli0 cations it suffices to consider the approximate distribution obtained by setting  = 1 (and hence ω = 1) in Eq. (19). This simplified distribution has the PGF of the form

where z a−1



m

−m + B(k, 1 + w−1 )z k w

G(z; m, w) = exp

∞

.

(20)

k=1

It is noteworthy that the beta function sequence B(k, 1 + w−1 ) in Eq. (20) can be computed with simple arithmetic operations. Specifically, one has B(1, 1 + w−1 ) = 1/(1 + w−1 ),

j−1 1 −1 1+w j + w−1 k

B(k, 1 + w−1 ) =

(k ≥ 2).

(21)

j=2

Eqs. (20) and (21) first appeared in Ref. [38]. One can compute the beta function sequence recursively: B1 = 1/(1 + w−1 ), Bk =

k−1 Bk−1 k + w−1

(22)

(k ≥ 2).

Therefore, the mutant probabilities are recursively calculated as follows: p(0; m, w) = e−m , m

jBj p(k − j; m, w) (k ≥ 1). kw k

p(k; m, w) =

(23)

j=1

The following h sequence defined in Ref. [45] is needed: h0 = −1, hk =

1 B(k, 1 + w−1 ) = w−1 Bk w

(k ≥ 1).

(24)

w=0.75

w=

w=1.00

Because both stains grow exponentially under laboratory conditions, Eq. (17) is an estimate of

5 × 10−9

1 × 10−8

1.5 × 10−8

2 × 10−8

Fig. 3. Ten thousand 15-tube fluctuation experiments were simulated with a mutation rate of 1.0 × 10−8 and relative fitness of 0.75. ML estimates of the mutation rate were generated by assuming w = 0.75 and w = 1.0 respectively. The distributions of the estimates of the mutation rate are summarized by boxplots. The vertical dashed line indicates the true mutation rate (1.0 × 10−8 ).

By means of sequence convolution, the first and second derivatives (with respect to m) of the mutant probabilities are computed as follows: p(1) (k; m, w) = hk ∗ p(k; m, w), p(2) (k; m, w) = hk ∗ p(1) (k; m, w).

(25)

Here the symbol * denotes the usual sequence convolution operator. The above discussion is a specialization of results given in Ref. [45]. ML estimates of m and likelihood-based CIs for m can be computed following the computational procedures described in Ref. [45] for the LC model. Likelihood ratio tests can be formed as outlined in Ref. [51]. Note that an alternative approach to accounting for fitness is to use information gleaned from a fitness assay to form a prior distribution of the parameter w. Bayesian inference about the parameter m can then be performed in a fashion similar to the procedure described in Ref. [48]. It has been generally believed that w < 1 for most mutant strains. As a corollary, most mutation rates were underestimated in the past 70 years or so, because w = 1 has been a tacit assumption in the analysis of fluctuation assay data. A more complex picture began to emerge recently. In the experiments of Scanlan et al. [36] w varies around 1.0. More remarkably, McElroy et al. [29] found that most mutations in their experiments caused large increases in w. Therefore, assessing the effects of mis-specification of w on the estimation of mutation rates is an important task. A simulation study was conducted in probe of this issue. Ten thousand 15-tube experiments were simulated using parameter values ˇmutant = 0.75, ˇwild-type = 1.0, N0 = 1, Nt = 3 ×108 and a mutation rate of 1.0 × 10−8 . ML estimates of the mutation rate were then generated by assuming w = 1.0 and 0.75 respectively. The distributions of these estimates are summarized in Fig. 3. The median of the estimated mutation rate under the assumption w = 0.75 is 1.01 × 10−8 , while under the assumption w = 1.0 it is 8.92 × 10−9 . The error due to failure to account for relative fitness is noticeable, although not alarmingly large. Another way to get a sense of the effects of mis-specification of w is the following semi-fictitious example. The experimental data of Demerec [11] were re-analyzed by assuming a sequence of values of w that were observed by Scanlan et al. [36]. ML estimates of the mutation rate along with 95% CIs are presented in Table 2. These values of w are likely to be encountered in other experiments. The estimates deviate from the benchmark value (obtained with w = 1) no more than 20%. This simple example

12

Q. Zheng / Mutation Research 781 (2015) 7–13

Table 2 Mutation rate estimate changes with fictitious fitness.

6. Discussion

Fitness w

Mutation rate

Lower limit

Upper limit

% Change

0.8435 0.8870 0.9971 1.0000 1.0496 1.0524 1.1175 1.1730 1.2173 1.2311 1.2357

6.52 6.28 5.72 5.71 5.48 5.47 5.19 4.97 4.81 4.76 4.74

5.30 5.07 4.57 4.55 4.35 4.34 4.09 3.89 3.75 3.71 3.69

7.83 7.56 6.96 6.94 6.70 6.68 6.38 6.14 5.96 5.90 5.88

14.31 10.04 0.24 0.00 3.98 4.19 9.05 12.90 15.80 16.67 16.96

The Demerec data [11] were re-analyzed by assuming a sequence of values of relative fitness as reported in Ref. [36]. Mutation rates (per 108 cell divisions) along with their 95% CIs were calculated using the algorithm discussed in Section 4. Percentage change (last column) refers to change in estimated mutation rate relative to that obtained by assuming w = 1.

suggests that, if w does not depart too much from unit, the LD model can yield useful estimates of mutation rates.

5. Comparison of mutation rates Comparison of mutation rates using the LD protocol is the most challenging aspect in the analysis of fluctuation assay data. To compare two mutation rates, two separate experiments are carried out. If average final cell population sizes in the two experiments, say Nt1 and Nt2 , are approximately the same, the comparison is straightforward, because comparing the two mutation rates is the same as comparing the two m parameters. The Mann–Whitney (MW) test is one obvious candidate statistical test. It has been shown [50] that likelihood ratio test (LRT) can be computed under this condition and is more powerful. When the condition Nt1 = Nt2 fails to hold, the MW test cannot be directly applied, but the empirical method of checking overlapping of CIs seems to still work [50]. Note that, even if the experimentalist can guarantee that the incubation periods of the two experiments are the same, the average final cell population sizes can differ appreciably due to differential growth rates of the two bacterium strains or some other factors. This difficult has been overcome recently [51]. A general framework for LRT has been developed that takes into consideration R, the ratio of Nt2 to Nt1 . Within this framework the two experiments can even differ in plating efficiency or relative fitness. With the LRT approach, multiple comparison can be performed, and the concept of false discovery rate [3] is directly applicable. Multiple comparison is an important technique to deal with a family of many experiments, e.g., the investigation reported in Ref. [41]. To show the usefulness of the new approach, consider the first experiment (strain H37Rv) and the last experiment (strain E 47/94) reported in Ref. [41]. If the MW test were applied directly to these two experiments, one would have p = 5.4 × 10−4 . However, this small p-value only suggests a difference between the expected numbers of mutations per culture in the two experiments, which may not be relevant to the comparison of the two mutation rates. In this case, m1 = 11.06 and m2 = 3.77, so the MW test rightfully yields a small p-value. However, when scaled by the final cell population sizes, the two mutation rates are 9.63 × 10−9 and 8.38 × 10−9 , respectively. Performing an LRT that accounts for the fact that R = 0.9 : 2.3 (and also the fact that  = 0.2 in both experiments), one obtains a p-value of 0.56. That is, the two mutation rates do not differ significantly, even though the expected numbers of mutations per culture (m) differ considerably. Improper use of the MW test can yield misleading results.

This paper aims toward updating and supplementing information given by a set of 5 papers that have served collectively as a practical guide to the use of the LD protocol. Important additions to the toolbox for studying microbial mutation rates are methods for adjusting for plating efficiency, for adjusting for relative fitness, and for comparison of mutation rates that can accommodate differences in final population size, plating efficiency and relative fitness. It also sounds a cautionary note about methods based on the sample median. The new methods advocated here were all derived using the likelihood principle. The likelihood approach can sometimes seem too trite to be exciting, but it remains the most reliable approach to tackling the three kinds of common problems in the analysis of fluctuation assay data—point estimation, interval estimation and comparison. The recent release of rSalvador [49], adapted from its predecessor SALVADOR [43], has made possible routine use by bench scientists of all the methods discussed here. The user only needs a minimal familiarity with the R language, a free computing platform favored increasingly by biologists. To begin with, the user can store the mutant count data from a fluctuation experiments in a single-column text file named, say, example1.txt. To make the data available to rSalvador, execute the R command y=import.text.data(‘example1.txt’) Now suppose that the text file contains the well-known experimental data of Demerec [11]. For the sake of illustration, assume that from a fitness assay the experimentalist learns that the relative fitness is w = 0.21, which is one of the fitness values observed by Billington et al. [4]. To get an ML estimate of m, execute the R command newton.MK(y,w=0.21) ˆ = 25.86 (under the assumption The user obtains an ML estimate m ˆ = 10.84). Furthermore, to set up a 95% CI for m, the w = 1.0, m experimentalist executes confint.MK(y,w=0.21) which yields a 95% CI of the form (22.98, 28.81). The difficulties discussed here, some appearing to be rather intricate, were caused by the experimentalist’s inability to observe mutations directly. If, when a wide-type cell divides, the genotypes of the two daughter cells could be observed and recorded, estimating a mutation rate in the laboratory would be all but trivial. Since the advent of the LD protocol, the barrier to actually seeing mutations occur has largely remained impenetrable. Encouraging advances have recently been made to break this barrier. For example, by fluorescently labeling mismatch repair proteins in Escherichia coli cells, Elez et al. [14] were able to count nascent base substitutions with the aid of time-lapse microscopy. This new technique for seeing mutations in living cells, while brilliant, is far from routinely applicable. The LD protocol remains the most practical approach for the estimation and comparison of microbial mutation rates. The present paper aims to help the experimentalist tap more fully the LD protocol’s potential by offering new or refined methods. Standardization of methods is not a practical remedy for all existing problems, because there are still open problems and because the field is still evolving. An alternative approach is that researchers voluntarily make their experimental data available to the reader. In the early days of fluctuation experiments, when the term “reproducible experiment” did not exist, experimental data were often published along with conclusions drawn from them.

Q. Zheng / Mutation Research 781 (2015) 7–13

In the past decade or so, this important scientific practice seemed to be on the verge of extinction. There is no denying that confusion about existing methods was partly responsible for some researchers’ reluctance to report raw experimental data. If raw experimental data were readily available (either presented in the paper or deposited as supplementary material), a diligent reader can easily verify major conclusions with the aid of a convenient software tool like rSalvador. In addition to the mutant count data, the researcher should also report N0 , Nt , and if applicable, plating efficiency, relative fitness, median cell life cycle and some relevant laboratory conditions like temperature. These data, like sequencing data, are of interest beyond the research project that led to the experiments. Acknowledgments I am indebted to three anonymous referees for detailed, constructive comments, which enabled me to restructure an earlier draft such that the resulting version is accessible to a wider audience. I am also profoundly grateful to G. Crouse, P. Foster and J. Jungck for encouraging me to launch the rSalvador project. Eric H. Zheng, then a computer engineering major at the University of Illinois at Urbana-Champaign, provided technical assistance in migrating rSalvador from Linux to the Windows platform. A major portion of the investigation relied on an IBM iDataPlex cluster and an IBM NeXtScale cluster, managed by Texas A&M University Supercomputing Facility. References [1] G. Asteris, S. Sarkar, Bayesian procedures for the estimation of mutation rates from fluctuation experiments, Genetics 142 (1996) 313–326. [2] H. Bachmann, M.J.C. Starrenburg, D. Molenaar, M. Kleerebezem, J.E.T. van Hylckama Vlieg, Microbial domestication signatures of Lactococcus lactis can be reproduced by experimental evolution, Genome Res. 22 (2012) 115–124. [3] Y. Benjamini, Y. Hochberg, Controlling the false discovery rate: a practical and powerful approach to multiple testing, J. R. Stat. Soc. Ser. B 57 (1995) 289–300. [4] O.J. Billington, T.D. MdHugh, S.H. Gillespie, Physiological cost of rifampin resistance induced in vitro in Mycobacterium tuberculosis, Antimicrob. Agents Chemother. 43 (1999) 1866–1869. [5] L. Boe, T. Tolker-Nielsen, K.-M. Eegholm, H. Spliid, A. Vrang, Fluctuation analysis of mutations to nalidixic acid resistance in Escherichia coli, J. Bacteriol. 176 (1994) 2781–2787. [6] J. Cairns, J. Overbaugh, S. Miller, The origin of mutants, Nature 335 (1988) 142–145. [7] F. Calcuttawala, C. Hariharan, G.P. Pazhani, S. Ghosh, T. Ramamurthy, Activity spectrum of colicins produced by Shigella sonnei and genetic mechanism of colicin resistance in conspecific S. sonnei strains and Escherichia coli, Antimicrob. Agents Chemother. 59 (2015) 152–158. [8] A. Couce, J. Blazquez, Estimating mutation rates in low-replication experiments, Mutat. Res. 714 (2011) 26–32. [9] A.J. Gruber, A.L. Erdem, G. Sabat, K. Karata, M.M. Jaszczur, D.D. Vo, T.M. Olsen, R. Woodgate, M.F. Goodman, M.M. Cox, A RecA protein surface required for activation of DNA polymerase V, PLoS Genet. 11 (3) (2015) e1005066. [10] J.A.G.M. de Visser, R.E. Lenski, Long-term experimental evolution in Escherichia coli. XI. Rejection of non-transitive interactions as cause of declining rate of adaptation, BMC Evol. Biol. 2 (2002) 19. [11] M. Demerec, Production of staphylococcus strains resistant to various concentrations of penicillin, Proc. Natl. Acad. Sci. U. S. A. 31 (1945) 16–24. [12] J.W. Drake, A constant rate of spontaneous mutation in DNA-based microbes, Proc. Natl. Acad. Sci. U. S. A. 88 (1991) 7160–7164. [13] S.F. Elena, R.W. Lenski, Evolution experiments with microorganisms: the dynamics and genetic bases of adaptation, Nat. Rev. Genet. 4 (2003) 457–469. [14] M. Elez, A.W. Murray, L.-J. Bi, X.-E. Zhang, I. Matic, M. Radman, Seeing mutations in living cells, Curr. Biol. 20 (2010) 1432–1437. [15] FALCOR, web page at http://www.keshavsingh.org/protocols/FALCOR.html. [16] R.C.E. Flohr, C.J. Blom, P.B. Rainey, H.J.E. Beaumont, Founder niche constrains evolutionary adaptive radiation, Proc. Natl. Acad. Sci. U. S. A. 110 (2013) 20663–20668. [17] P.L. Foster, Methods for determining spontaneous mutation rates, Methods Enzymol. 409 (2006) 195–213. [18] B.M. Hall, C.-X. Ma, P. Liang, K.K. Singh, Fluctuation AnaLysis CalculatOR: a web tool for the determination of mutation rate using Luria–Delbrück fluctuation analysis, Bioinformatics 25 (2009) 1564–1565.

13

[19] M.E. Jones, S.M. Thomas, A. Rogers, Luria–Delbrück fluctuation experiments: design and analysis, Genetics 136 (1994) 1209–1216. [20] W.S. Kendal, P. Frost, Pitfalls and practice of Luria–Delbrück fluctuation analysis: a review, Cancer Res. 48 (1988) 1060–1065. [21] A.L. Koch, Mutation and growth rates from Luria–Delbrück fluctuation tests, Mutat. Res. 95 (1982) 129–143. [22] R. Kraˇsovec, R.V. Belavkin, J.A.D. Aston, A. Channon, E. Aston, B.M. Rash, M. Kadirvel, S. Forbes, C.G. Knight, Mutation rate plasticity in rifampicin resistance depends on Escherichia coli cell–cell interactions, Nat. Commun. 5 (2014) 3742, http://dx.doi.org/10.1028/ncomms4747. [23] E.A. Lea, C.A. Coulson, The distribution of the numbers of mutants in bacterial populations, J. Genet. 49 (1949) 264–285. [24] D.L. Lindstrom, D.E. Gottschling, The mother enrichment program: a genetic system for facile replicative life span analysis in Saccharomyces cerevisiae, Genetics 183 (2009) 413–422. [25] S.E. Luria, M. Delbrück, Mutations of bacteria from virus sensitivity to virus resistance, Genetics 28 (1943) 491–511. [26] S.E. Luria, A Slot Machine, A Broken Test Tube: An Autobiography, Harper & Row, New York, 1984. [27] W.T. Ma, G.v.H. Sandri, S. Sarkar, Analysis of the Luria and Delbrück distribution using discrete convolution powers, J. Appl. Prob. 29 (1992) 255–267. [28] B. Mandelbrot, A population birth-and-mutation process, I: Explicit distributions for the number of mutants in an old culture of bacteria, J. Appl. Prob. 11 (1974) 437–444. [29] K.E. McElroy, J.G.K. Hui, J.K.K. Woo, A.W.S. Luk, J.S. Webb, S. Kjelleberg, S.A. Rice, T. Thomas, Strain-specific parallel evolution drives short-term diversification during Psedomonas aeruginosa biofilm formation, Proc. Natl. Acad. Sci. U. S. A. 111 (2014) E1419–E1427. [30] H.B. Newcombe, Delayed phenotypic expression of spontaneous mutations in Escherichia coli, Genetics 33 (1948) 447–476. [31] S. Oide, W. Gunji, Y. Moteki, S. Yamamoto, M. Suda, T. Jojima, H. Yukawa, M. Inui, Thermal and solvent stress cross-tolerance conferred to Corynebacterium glutamicum by adaptive laboratory evolution, Appl. Environ. Microbiol. 81 (2015) 2284–2298. [32] H.A. Orr, Fitness and its role in evolutionary genetics, Nat. Rev. Genet. 10 (2009) 531–539. [33] C.F. Pope, D.M. O’Sullivan, T.D. McHugh, S.H. Gillespie, A practical guide to measuring mutation rates in antibiotic resistance, Antimicrob. Agents Chemother. 52 (2008) 1209–1214. [34] C. Rajanna, G. Quellette, M. Rashid, A. Zemla, M. Karavis, C. Zhou, T. Revazishvili, B. Redmond, L. McNew, L. Bakanidze, P. Imnadze, B. Rivers, E.W. Skowronki, H.P. O’Connell, A. Sulakvelidze, H.S. Gibbons, A strain of Yersinia pestis with a mutator phenotype from the Republic of Georgia, FEMS Microbiol. Lett. 343 (2013) 113–120. [35] W.A. Rosche, P.L. Foster, Determining mutation rates in bacterial populations, Methods 20 (2000) 4–17. [36] P.D. Scanlan, A.R. Hall, G. Blackshields, V.-P. Friman, M.R. Davis Jr., J.B. Goldberg, A. Buckling, Coevolution with bacteriophage drives genome-wide host evolution and constrains the acquisition of abiotic-beneficial mutations, Mol. Biol. Evol. 32 (2015) 1425–1435. [37] A. Srivastava, D. Degen, Y.W. Ebright, R.H. Ebright, Frequency, spectrum, and nonzero fitness costs of resistance to myxopyronin in Staphylococcus aureus, Antimicrob. Agents Chemother. 56 (2012) 6250–6266. [38] F.M. Stewart, D.M. Gordon, B.R. Levin, Fluctuation analysis: the probability distribution of the number of mutants under different conditions, Genetics 124 (1990) 175–185. [39] F.M. Stewart, Fluctuation tests: how reliable are the estimates of mutation rates, Genetics 137 (1994) 1139–1146. [40] D.M. Stoebel, C.J. Dorman, The effect of mobile element IS10 on experimental regulatory evolution in Escherichia coli, Mol. Biol. Evol. 27 (2010) 2105–2112. [41] J. Werngren, S.E. Hoffner, Drug-susceptible Mycobacterium tuberculosis Beijing genotype does not develop mutation-conferred resistance to rifampin at an elevated rate, J. Clin. Microbiol. 41 (2003) 1520–1524. [42] Q. Zheng, Progress of a half century in the study of the Luria–Delbrück distribution, Math. Biosci. 162 (1999) 1–32. [43] Q. Zheng, Statistical and algorithmic methods for fluctuation analysis with SALVADOR as an implementation, Math. Biosci. 176 (2002) 237–252. [44] Q. Zheng, Update on estimation of mutation rates using data from fluctuation experiments, Genetics 171 (2005) 861–864. [45] Q. Zheng, New algorithms for Luria–Delbrück fluctuation analysis, Math. Biosci. 196 (2005) 198–214. [46] Q. Zheng, On Haldane’s formulation of Luria and Delbrück’s mutation model, Math. Biosci. 209 (2007) 500–513. [47] Q. Zheng, A note on plating efficiency in fluctuation experiments, Math. Biosci. 216 (2008) 150–153. [48] Q. Zheng, A. Bayesian approach for correcting for partial plating in fluctuation experiments, Genet. Res. 93 (2011) 351–356. [49] Q. Zheng, rSalvador: an R tool for the Luria–Delbrück fluctuation assay, 2014 http://eeeeeric.github.io/rSalvador. [50] Q. Zheng, Methods for comparing mutation rates using fluctuation assay data, Mutat. Res. 777 (2015) 20–22. [51] Q. Zheng, Comparing mutation rates under the Luria–Delbrück protocol, 2015 (in review).