Our research groups have previously published a study estimating the bottleneck size of Cauliflower mosaic virus (CaMV) during host plant colonization [1]. Two methods were used in that study, both based on the temporal evolution of neutral marker frequencies; one uses estimates of FST, and the other tracks changes in marker frequency variance over time. Here we report that (i) the variance-based method was actually published by Felsenstein [2], a reference we were not aware of when publishing the original study; (ii) equation (4) in [1] is actually an approximation; and (iii) these methods rely on the assumption that the bottleneck size is constant, which can be relaxed by assuming it is a Poisson-distributed random variable. Equation (4) in [1] reads Nv = p(1−p)/[Var(p’)−Var(p)], where Nv is the estimated bottleneck size, and p and p’ represent the frequency of a neutral marker before and after the bottleneck, respectively. The approximation in this equation concerns the numerator. Its exact formulation should be E(p(1 − p)) and not E(p)(1 − E(p)), which is the expression implied in [1] and used in the numerical application on CaMV (E denotes the expected value). Based on the definition of the variance, it can be shown that E(p(1 − p)) = E(p)(1 − E(p)) − Var(p), which is equivalent to the numerator of equation (14) in [2]. Thus, the approximation will give relatively accurate results whenever Var(p) is negligible relative to E(p)(1 − E(p)). This was indeed the case in the study by Monsion et al. [1], as can be seen in Table 1, which compares the Nv values published in [1] to those obtained when applying the exact expression. The results of that paper are thus not qualitatively affected and, consequently, its conclusions remain unchanged. This method was also used in another study [3], the results of which are also marginally affected quantitatively, but not at all qualitatively (Table 2). To illustrate the application of these methods in a situation in which the estimated bottleneck size is lower, we used data from Tromas et al. (Table 3) [4]. Table 1 Comparison of estimates of Nv using the approximate formula published in [1] (values in Table 2 of this publication), those using the exact expression, and those assuming that bottleneck size is a zero-truncated-Poisson-distributed random variable. Table 2 Comparison of estimates of Nv published in [3] using the approximate formula (values in table S1 of [3]), those using the exact expression, and those assuming that bottleneck size is a zero-truncated-Poisson-distributed random variable. Table 3 Estimates of Nv using data published in [4] for leaf level 5. Both Felsenstein’s formula [2] and its approximation assume that the bottleneck size N is constant (i.e., identical in all the sampled plants). However, we can make the more realistic assumption that N follows a zero-truncated Poisson (ZTP) distribution (i.e., N can vary among experimental replicates according to a Poisson distribution, but p’ cannot be measured when N = 0 and the plant is discarded, hence the zero-truncation): N∼ZTP(nvP),so∀k∈N*,P(N=k)=nvPkk!(envP−1). Among the N genomes that go through the bottleneck, the number X bearing the neutral marker is drawn according to its pre-bottleneck frequency p from the binomial distribution: X~B(N,p). Thus, the marker frequency after the bottleneck is p’ = X/N. The variance of this marker frequency can be expressed as: Var(p’) = Var[E(p’|N,p)]+E[Var(p’|N,p)]. Because E(p’|N,p) = E(X|N,p)/N and Var(p’|N,p) = Var(X|N,p)/N2, and because N and p are independent, we then get: Var(p’) = Var(p)+E[p(1−p)]×E(1/N). Using the properties E(1/N)=∑k=1∞(1kP(N=k)) and ∑k=1∞nvPkk!k=Ei(nvP)−ln(nvP)−γ, where Ei=∫−∞xettdt is the exponential integral function and γ the Euler-Mascheroni constant, one can estimate the bottleneck size nvP by numerically solving the following equation: E(p)(1−E(p))−Var(p)Var(p')−Var(p)=envP−1Ei(nvP)−ln(nvP)−γ. We recommend using this estimation procedure when Nv values are below 20, i.e., when the relative estimation error is higher than 5%. The following lines of R code may be used to estimate nvP: library(gsl) EulerMaschCst<−(−digamma(1)) p<−c(…,.,…) # vector of initial frequencies p_prime<−c(…,.,…) # vector of final frequencies funestim<−function(n_vP) (mean(p)*(1−mean(p))−var(p))/(var(p_prime)−var(p)) − (exp(n_vP)−1)/(expint_Ei(n_vP)−log(n_vP)−EulerMaschCst) uniroot(funestim,interval = c(1e−6,600))
Read full abstract
May 12, 2021
Jharna Pradhan + 1
Open Access