A large effective population size for established withinhost influenza virus infection
Abstract
Strains of the influenza virus form coherent global populations, yet exist at the level of single infections in individual hosts. The relationship between these scales is a critical topic for understanding viral evolution. Here we investigate the withinhost relationship between selection and the stochastic effects of genetic drift, estimating an effective population size of infection N_{e} for influenza infection. Examining wholegenome sequence data describing a chronic case of influenza B in a severely immunocompromised child we infer an N_{e} of 2.5 × 10^{7} (95% confidence range 1.0 × 10^{7} to 9.0 × 10^{7}) suggesting that genetic drift is of minimal importance during an established influenza infection. Our result, supported by data from influenza A infection, suggests that positive selection during withinhost infection is primarily limited by the typically short period of infection. Atypically long infections may have a disproportionate influence upon global patterns of viral evolution.
Introduction
The evolution of the influenza virus may be considered across a broad range of scales. On a global level, populations exhibit coherent behaviour (Buonagurio et al., 1986; Fitch et al., 1997; Bedford et al., 2015), evolving rapidly under collective host immune pressure (Ferguson et al., 2003; Grenfell et al., 2004). On another level, these global populations are nothing more than very large numbers of individual host infections, separated by transmission events.
Despite the clear role for selection in global influenza populations, recent studies of withinhost infection have suggested that positive selection does not strongly influence evolution at this smaller scale (Debbink et al., 2017; McCrone et al., 2018; Han et al., 2019). Contrasting explanations have been given for this, with suggestions either that selection at the withinhost level is intrinsically inefficient, being dominated by stochastic processes (McCrone et al., 2018), or that while selection is efficient, a mismatch in timing between the peak viral titre and the host adaptive immune response prevents selection from taking effect (Han et al., 2019).
To resolve this issue, we evaluated the relative importance of selection and genetic drift during a case of influenza infection. The balance between these factors is determined by the effective size of the population, denoted N_{e}. If N_{e} is high, selection will outweigh genetic drift, even where differences in viral fitness are small (Rouzine et al., 2001). By contrast, if N_{e} is low, less fit viruses are more likely to outcompete their fitter compatriots.
Estimating N_{e} is a difficult task, with a long history of method development in this area (Wright, 1938; Wang et al., 2016; Khatri and Burt, 2019). A simple measure of N_{e} may be calculated by matching the genetic change in allele frequencies in a population with the changes occurring in an idealised population evolving under genetic drift (Kimura and Crow, 1963). However, such estimates are vulnerable to distortion, for example being reduced by the effect of positive selection in a population. Where the global influenza A/H3N2 population is driven by repeated selective sweeps (Fitch et al., 1991; Rambaut et al., 2008; Strelkowa and Lässig, 2012) a neutral estimation method suggests a value for N_{e} not much greater than 100 (Bedford et al., 2010). While methods for jointly estimating N_{e} and selection exist, they are limited in considering only a few loci in linkage disequililbrium (Bollback et al., 2008; Feder et al., 2014; Foll et al., 2014; Terhorst et al., 2015; Rousseau et al., 2017). Nontrivial population structure can affect estimates (Laporte and Charlesworth, 2002); a growing body of evidence supports the existence of structure in withinhost influenza infection (Lakdawala et al., 2015; Sobel Leonard et al., 2017a; Richard et al., 2018; Hamada et al., 2012). While careful experimental techniques can reduce sequencing error (McCrone and Lauring, 2016), noise from sequencing and unrepresentative sample collection combine (Illingworth et al., 2017), potentially confounding estimates of N_{e} in viral populations (Lumby et al., 2018). If N_{e} is high, any signal of drift can be obscured by noise.
We here estimate a mean effective population size for an established withinhost influenza B infection using data collected from a severely immunocompromised host. While the viral load of the infection was not unusual for a hospitalised childhood infection (Wishaupt et al., 2017), an absence of cellmediated immunity led to the persistence of the infection for several months (Lumby et al., 2020). Given extensive sequence data collected during infection, the reduced role of positive selection, combined with novel methods to account for noise and population structure, enabled an improved inference of N_{e}. The large effective size we infer suggests that selection acts in an efficient manner during an established influenza infection. Even in more typical cases, the influence of positive selection is likely to be limited only by the duration of infection.
Results and discussion
Viral samples were collected at 41 time points spanning 8 months during the course of an influenza B infection in a severely immunocompromised host (Figure 1A). Clinical details of the case, and the use of viral sequence data in evaluating the effectiveness of clinical intervention, have been described elsewhere (Lumby et al., 2020). After unsuccessful treatment with oseltamivir, zanamivir and nitazoxanide, a bone marrow transplant and favipiravir combination therapy led to the apparent clearance of infection. Apart from a single exception, biweekly samples tested negative for influenza across a period of close to two months. A subsequent resurgence of zanamivirresistant infection was cleared by favipiravir and zanamivir in combination.
Phylogenetic analysis of wholegenome viral consensus sequences showed the existence of nontrivial population structure, with at least two distinct clades emerging over time (Figure 1B, Figure 1—figure supplement 1); we term these clades A and B. Having diverged, the two clades persisted across several months of infection. Haplotype reconstruction showed that samples from clade B were comprised of distinct viral haplotypes to those from clade A; similar patterns were observed in different viral segments (Figure 1—figure supplement 2). The October 4^{th} sample is intermediate between the initial and final samples collected (Figure 1D). We suggest that, from a common evolutionary origin, Clade B slowly evolved away from the initial consensus, while viruses in clade A stayed close in sequence space to this consensus. The cladal structure suggests the existence of spatially distinct viral populations in the host, samples stochastically representing one population or the other.
To estimate the effective population size, we analysed genomewide sequence data from samples in clade A collected before first use of favipiravir. A method of linear regression was used to quantify the rate of viral evolution, measuring the genetic distance between samples as a function of increasing time between dates of sample collection. We inferred a rate equivalent to 0.051 substitutions per day (97.5% confidence interval 0.034 to 0.068) (Figure 2A), equivalent to 7.94 substitutions genomewide across 157 days of evolution. The vertical intercept of this line provides an estimate of the contribution of noise to the measured distance between samples, potentially arising from sequencing error or undiagnosed population structure. The identified value of close to 40 substitutions is equivalent to a betweensample allele frequency difference of approximately +/ 0.3% per locus. While considerable noise affects each sample, the dataset as a whole provides a clear signal of evolutionary change.
A simulation based analysis, measuring the extent of evolution in idealised WrightFisher populations (Kimura and Crow, 1963), inferred an effective population size of 2.5 × 10^{7} (95% confidence range 1.0 × 10^{7} to 9.0 × 10^{7}) for viruses in clade A before the use of favipiravir (Figure 2B). This value is substantially larger than estimates made recently for withinhost HIV infection (Pennings et al., 2014; Rouzine et al., 2014), and suggests that even weak selection could easily overcome genetic drift. Data from clade B gave a lower estimated value of 2 × 10^{6}, (95% confidence range 4 × 10^{5} to 2 × 10^{8}) perhaps reflecting the less frequent observation of samples in that clade (Figure 2C,D), and the bottleneck induced by favipiravir, which was spanned by the data used in this calculation.
Our value of N_{e} is representative of the population after the initial establishment of infection; the initial expansion of the viral population was not represented in our data. Population structure during the infection might have lowered the value we obtain (Whitlock and Barton, 1997). The partial onset of zanamivir resistant alleles (Jackson et al., 2005), sporadically observed at intermediate frequency in clade A after the administration of the drug (Figure 2—figure supplement 1), is suggestive of sampling a random mixture of viruses from resistant and susceptible subpopulations.
Our method equates change in a population with genetic drift (Kimura and Crow, 1963), neglecting the role of selection. As such, the influence of positive selection might have led us to underestimate N_{e}. While viral evolution was generally not driven by selection (Figure 2—figure supplement 2), positive selection (e.g. for zanamivir resistance) would increase the rate of viral evolution, lowering our inferred value. Selection may have influenced the division between clades, perhaps through the adaptation of the virus to specific local environments. Purifying selection may also have influenced the population in ways not accounted for by our method. Yet our result is clear. Once an infection is established, selection will dominate the stochastic effects of drift upon withinhost evolution.
The dataset we considered is particularly suited to our calculation. The long period of infection combined with frequent sampling allowed for the characterisation of a slow rate of evolution amidst population structure and noise in the data. Further, the absence of strong selection reduced the error in our inference approach, which assumed an idealised neutral population. To provide further validation we repeated our approach on data describing longterm influenza A/H3N2 infection in four immunocompromised adults (Xue et al., 2017). The estimates for N_{e} we obtained, of between 3 × 10^{5} and 1 × 10^{6} (Figure 2—figure supplement 3), while high, were smaller than for our flu B case, potentially being reduced by an increased influence of selection.
We believe that our study provides a first realistic estimate of withinhost effective population size for severe influenza infection in humans. The viral load in the influenza B case was high, representative of hospitalised cases of childhood influenza infection. However, the magnitude of our inferred effective size, of order 10^{7}, suggests that selection will predominate over drift even in more typical cases. Mean CT values for influenza in nonhospitalised children have been reported as around 10 units lower than those for hospitalised cases (Wishaupt et al., 2017); an order of magnitude calculation suggests an Ne, upon the establishment of infection, of approximately 10^{4} in such cases. Such a value again reflects an established population, not accounting for the initial population bottleneck. It has the implication that the evolution of a measurable variant (i.e. at a frequency of 1% or above) will be dominated by selection of a magnitude of 1% or greater per generation (Rouzine et al., 2001).
Our result supports the idea that a tight transmission bottleneck (McCrone et al., 2018; Valesano, 2020; Ghafari et al., 2020) followed by a short period of infection is sufficient to explain the observed lack of withinhost variation in typical cases of influenza (Debbink et al., 2017; McCrone et al., 2018); the stochastic effects of genetic drift do not limit the impact of positive selection. Variants arising through de novo mutation would require strong selection to reach a substantial frequency during infection (Zhao et al., 2019), particularly if the onset of selection is delayed (Miao et al., 2010; Illingworth et al., 2014; Morris, 2020). We suggest that, while not being confounded by drift, selection does not usually have time to fix novel variants in the population, exceptions including the emergence of antiviral resistance and some cases of longer infection (Xue et al., 2017; Gubareva et al., 1998; Snydman, 2006; Centers for Disease Control and Prevention (CDC), 2009; Imai et al., 2020; Rogers et al., 2015).
Our result highlights the potential importance of longer infections in the adaptation of global influenza populations, particularly where some adaptive immune response remains. A newly emergent variant under strong positive selection increases faster than linearly in frequency (Haldane, 1924). Given a large N_{e}, implying efficient selection, additional days of infection will have a disproportionate influence upon the potential transmission of adaptive variants. This does not imply that longer infections are the sole driving force behind global viral adaptation; selective effects affecting viral transmissibility (Lumby et al., 2018) would provide an alternative explanation. However, our work suggests that longerterm infections may be an important area of study in the quest to better understand global influenza virus evolution.
Materials and methods
Summary
In a singlelocus haploid system, the expected change in a variant allele with frequency q caused by genetic drift is given by the formula (Charlesworth, 2009)
This fact has been exploited to evaluate the size of transmission bottlenecks in influenza infection, comparing statistics of genome sequence data collected before and after a transmission event (Poon et al., 2016; Sobel Leonard et al., 2017b). Such a calculation may be affected by noise in the sampling or sequencing of a population, particularly where the extent of noise outweighs the genuine change in a population (Lumby et al., 2018). Here we suggest that, given multiple samples from a population, an alternative approach is possible; we use this to derive a more robust estimate of N_{e}. By means of evolutionary simulations we estimate N_{e} for cases of withinhost influenza infection.
Sequence data and bioinformatics
Request a detailed protocolSequence data describing the evolution of the infection was generated as part of a previous study (Lumby et al., 2020). Data, edited to remove human genome sequence data, have been deposited in the Sequence Read Archive with BioProject ID PRJNA601176. The HCV data used in validating the sequencing pipeline (see below) were previously deposited in the Sequence Read Archive with BioProject ID PRJNA380188. Processed files describing raw variant frequencies for both datasets are available, along with code used in this project, at https://github.com/cjri/FluBData (copy archived at https://github.com/elifesciencespublications/FluBData; Illingworth, 2020a).
Shortread data were aligned first to a broad set of influenza sequences. Sequences from this set to which the highest number of reads aligned were identified and used to carry out a second shortread alignment. The SAMFIRE software package was then used to filter the shortread data with a PHRED score cutoff of 30, to identify consensus sequences, and to calculate the number of each nucleotide found at each position in the genome. SAMFIRE is available from https://github.com/cjri/samfire (Illingworth, 2020b).
Calculation of evolutionary distances
Request a detailed protocolVariant frequencies at different time points during infection were used to calculate a rate of change in the population over time. We define q(t) as a 4 x L element vector describing the frequencies of each of the nucleotides A, C, G, and T at each locus in the viral genome at time t. We next define a distance between vectors q. Considering a single locus in the genome, we calculate the change in allele frequencies over time via a generalisation of the Hamming distance
where the term inside the sum indicates the absolute difference between the frequency of allele a at locus i. The statistic d_{i} is equal to one in the case of a substitution, for example where only A nucleotides are observed in one sample and only G nucleotides in another. However, in contrast to the Hamming distance it further captures smaller changes in allele frequencies, lesser changes producing values between zero and one, such that a change of a variant frequency from 45% to 55% at a twoallele locus would equate to a distance of 0.1, representing half of the sum of the absolute changes in each of the two frequencies. The total distance between the two vectors may now be calculated as
where the sum over i is conducted over all loci in the viral genome.
Sequence distances for nonsynonymous and synonymous mutations were calculated in a similar manner, with the exception that distances were calculated over individual nucleotides rather than in a perlocus manner. We calculated
and
where A_{N,i} and A_{S,i} are the sets of nucleotides a and positions i in the genome which respectively induce nonsynonymous and synonymous changes in the consensus sequence. Synonymous and nonsynonymous variants were identified with respect to influenza B protein sequences; a nucleotide substitution was defined as being nonsynonymous if it induced a change in the coded protein in at least one viral protein sequence. By contrast to our primary distance measurement, values for synonymous and nonsynonymous sites were calculated as mean distances per nucleotide, reflecting the differing numbers of each type of potential substitution in the viral genome.
Estimation of effective population size
Request a detailed protocolWe converted our measurements of sequence distance into an estimate of N_{e} by means of a simplified evolutionary model, assuming that all of the change in the population results from genetic drift. We first note the effect of error in measurements of the population upon our distance metric.
We suppose that at the time t, we make the observation:
where e is the error in measuring the population. Our definition of ‘error’ here is a broad one; we include both the potential for viral material in a single swab to not fully capture the entire viral diversity within the host and the potential for the sequencing pipeline to distort the composition of the material in the swab (Illingworth et al., 2017). In our distance calculation, we now have:
where the terms e_{i} are locusspecific errors in the measurement of allele frequencies; we write this equation in the form:
where E is the deviation incurred from the true distance.
Here, given only two errorprone samples from a system, separation of the real population distance and the error term is impossible. However, given multiple samples, an approximate separation can be made. We here use linear regression to fit a model to the observed distances, fitting the model:
for constant values k, approximating the rate of evolutionary change in the population per unit time, and E, approximating the mean amount of error in a measurement; here the term in vertical brackets is the absolute difference in time between samples i and j. This approach makes two approximations, which we believe to be either reasonable or possible to account for. Firstly, the model assumes that a linear model is appropriate to describe the change in the population over time; within our drift framework this is correct if the effective population size N_{e} is constant, and if the distribution of allele frequencies does not change over time. In our data, the consensus population declines approximately eightfold (Lumby et al., 2020), then undergoes a bottleneck due to the influence of favipiravir; we infer a representative mean value of N, selecting for clade A only samples collected before the bottleneck. Secondly, our model assumes that the deviation from truth in our distance metric does not change in a manner that is systematically associated with the time between samples. Regarding the sequencing process we believe this to be correct in so far as a consistent sequencing pipeline was used throughout. Regarding withinhost population structure we note in our data a divergence over time between samples from clade A and clade B, but split these samples to obtain distinct estimates of N_{e} for each clade. We note that large deviations from our model assumptions can be qualitatively identified by a poor fit between a simple regression model and the data.
Linear regression was performed using the Mathematica 11 software package, using the same package to calculate a 97.5% confidence interval for the calculated gradient, k.
WrightFisher simulation
Request a detailed protocolWe next approximated the behaviour of our system using a WrightFisher model, rewriting the first component of Equation 9 as
Here ΔD is a stochastic function describing the change in the population, measured according to the metric D, that arises from a single generation of genetic drift in a population with effective size N_{e} and initial allele frequencies q(t_{1}). Regarding these allele frequencies we note that the distribution of minor allele frequencies across the genome was reasonably constant between samples for which a good read depth was achieved (Figure 2—figure supplement 4; read depths for these data have previously been reported Lumby et al., 2020). To account for variance in these statistics we used different samples to initiate our simulations, reporting error bars across choices of q(t_{1}).
Our WrightFisher model simulated the evolution of the viral population for a single generation. Rates of evolution calculated from the sequence data were rates of change per day whereas a WrightFisher simulation gives an estimated rate of evolution per generation. We therefore scaled the former to match the experimentally ascertained estimate of 10 hr per generation for influenza B (Nobusawa and Sato, 2006).
To conduct a simulation we constructed a population of N viruses. Each simulated virus had a genome comprised of eight segments, each identical in length to the corresponding segment of the influenza B virus sampled from the patient. Observations from the clinical viral population were used to specify the genetic composition of the viral population at the beginning of the simulation. A simulated population of viral genomes was established. For each viral segment, a clinical sample was chosen at random. Nucleotide frequencies at each locus in the clinical sample (modified as described below) were used to generate a multinomial sample of viruses from the simulated population, assigning alleles to viruses in the simulated population according to the random sample. This step was repeated for each locus in the segment, with no intrinsic association between alleles at different loci. The sample collected on 30th November 2017 was excluded as a starting point from this analysis due to its low read depth; all other samples had a mean read depth in excess of 2000fold coverage.
Simulation of the population was conducted at the genomewide level. We simulated a single generation of the evolution of our population under genetic drift, generating a random sample of N whole viral genomes from the population. Intrasegment recombination was assumed to be negligible (Boni et al., 2008), while reassortment between segments was neglected in line with evidence from cases of human infection (Sobel Leonard et al., 2017a). We collected allele frequency data from the initial and final populations, using these to calculate the distance in sequence space through which the population had evolved according to the modified Hamming distance described above.
For each population size tested, our simulation was run 400 times, using the data to produce a 97.5% confidence interval for the extent of evolutionary change at a given effective population size. For each of these 400 replicate simulations, an independent random set of samples was chosen to initiate each of the eight simulated viral segments. The extent of evolution of the real population was compared to the results from our simulated populations, giving an inference of the effective size of the viral population.
Amendments were made to the above approach.
Accounting for falsepositive variants in sequencing: Estimating a false positive rate
Request a detailed protocolThe evolutionary distance ΔD(N,q(t_{1})) calculated by our method is dependent upon the vector of allele frequencies q. Given a greater number of polymorphic alleles in a system, the evolutionary distance, calculated as the sum of allele frequency changes, will also increase. While the experimental pipeline we used has been shown to perform well in capturing withinhost viral diversity (STOPHCV Consortium et al., 2016), the possibility remains that sequencing could contribute additional diversity to the initial populations used in our simulation. We therefore made an estimate of the extent to which our sequencing process led to the false identification of variants. To achieve this, we used data from a previous study describing the repeat sequencing of hepatitis C virus (HCV) samples from a host (Illingworth et al., 2017); data in this previous study were collected using the same sequencing pipeline as that used to collect the data considered here and therefore provide a generic measure of the level of false positive variation. The data we analysed, coded as HCV01 in the original study, comprised four clinical HCV samples, each of which was split following nucleic acid extraction. Some replicate samples were processed using a DNase depletion method before all samples went through cDNA synthesis, library preparation and sequencing. DNase depletion led to samples with lower read depth; we here compared sequence data collected from the nondepleted replicates of each sample. Variant frequencies within this dataset, where variation was observed in more than one sample, are shown in Figure 2—figure supplement 5.
Considering the real viral sample, we note that at any given genetic locus, a minority variant either exists or does not exist according to some welldefined criterion. (For the moment the way in which variation is defined is not important; methods for defining variation, which include the use of a frequency threshold, are discussed later.) We denote the possible states of a locus as P and N, according to whether the locus is positive or negative for variation. We suppose that the probability that a random locus in the genome has a minority variant is given by P_{P}, leading to the equivalent statistic P_{N} = 1 P_{P}.
Sequencing of a specific position in the genome results in the observation or nonobservation of a variant. In our data we have sets of two replicate observations of each position in the genome, giving for each minority variant the possible outcomes VV, VX, XV, and XX, where V corresponds to the observation of a variant, and X corresponds to the nonobservation of a variant. These observations contain errors; we denote the true positive, false positive, true negative and false negative rates of the variant identification process by P_{VP}, P_{VN}, P_{XN}, and P_{XP} respectively. In this notation, VP indicates the observation of a variant conditional on the variant being a true positive.
The underlying purpose of our calculation is to remove falsely detected variation from the population. We begin by assuming that the false negative rate of detecting variants is equal to zero. That is, where we do not see a variant in the sequence data, we assume that a variant is never actually present. This is a conservative step in so far as we never add unobserved variation to the population. Our assumption gives the result that the false negative rate, P_{XP} = 0. In so far that a variant is never unobserved it follows that the true positive rate P_{VP} = 1.
We may now construct expressions for the probabilities of observing each of the four possible outcomes. Noting that P_{VN} + P_{XN} = 1 we obtain
Thus the outcome probabilities may be expressed in terms of the underlying probability of a position having a variant, P_{P}, and the false positive rate P_{VN}.
We next processed our sequence replicate data, considering only sites that were sequenced to a read depth of at least 2000fold coverage. For each locus in a dataset, we calculated the observed frequency of each of the nucleotides A, C, G, and T, generating pairs which described these frequencies in each of our two replicate datasets. Removing pairs in which an allele has a frequency of more than 0.5 in either of the two datasets, we obtained a list of minority variants from each locus, generally comprising three allele frequency pairs per locus. If it is correct that two of the three minority alleles have very low frequencies, the frequencies are close to being statistically independent; the existence of a very few alleles of one minority type does not greatly affect the probability of another variant allele being observed in another read. We note that, of the more than 73 thousand sites sequenced, only 56, fewer than 0.1%, had more than one minority variant at a frequency greater than 1%. We proceeded on the assumption that each pair of minority frequencies was statistically independent of the others.
From the repeated observations of sites, we may count the number of observations of each of the four outcomes; given a total of N pairs we denote these as N_{VV}, N_{VX}, N_{XV}, and N_{XX}. Under our model of independent pairs we constructed the multinomial log likelihood of the underlying variant and false positive rates.
where the terms P_{ab} are constructed from P_{P} and P_{VN} according to the equations above.
Given a set of paired observations, we calculated the maximum likelihood values of P_{P} and P_{VN}. From these statistics we are able to calculate the positive predictive value of sequencing, namely the proportion of observed variants that are true positives. This is achieved by dividing the probability that a true positive was detected (equal to the number of true positives as P_{VP} = 1), by the probability that a variant was detected:
Frequency dependence of falsepositive variant calling
Request a detailed protocolWithin our data, our expectation was that minority variants at higher allele frequencies would be more likely to be observed as variants in both replicate samples. We note that, where a frequency cutoff is applied to identify variants, care is required in the above protocol. For example, if a hard threshold was applied, in which variants were called at 1% frequency, a variant that was detected at frequencies of 1.01% and 0.99% would be regarded as having been observed in one case, and not observed in the other, although it likely represents a consistent observation.
In order to assess the frequency dependence of our true positive rate, we defined minimum and maximum variant frequency thresholds q^{min} and q^{max}, and denoted the replicate observations of a minority variant frequency as q^{A} and q^{B} in the two samples. We further defined the frequency q^{cut} according to the formula:
We then defined regions of frequency space as follows:
These inequalities are illustrated in Figure 2—figure supplement 6.
In the above, q^{cut} functions to slightly harshen the criteria for detecting variants at low frequencies. If a variant is observed in one sample at frequency greater than q^{min}, then if q^{min} is greater than 0.2%, the frequency in the second sample had to be at least half q^{min} to be counted. If q^{min} was between 0.1% and 0.2%, the frequency in the second sample had to be at least 0.1%, while if q^{min} was less than 0.1%, the frequency in the second sample had to be at least q^{min}.
For different ranges of frequency values, q^{min} and q^{max}, the proportion of observed variants that were true positives was calculated according to the maximum likelihood method above, using these categorisations. Results are shown in Figure 2—figure supplement 7. In the process of setting up the initial state of our WrightFisher simulation variants observed in the sequence data were considered in turn, drawing a Bernoulli random variable for each variant. Variants were included in the initial simulated population with probability equal to the proportion of observed variants that were estimated to be true positives.
Accounting for mutationselection balance
Request a detailed protocolTo account for our neglect of mutation, a frequency cutoff was applied to our simulation data. Under a pure process of genetic drift, lowfrequency variants in our population are likely to die out, reaching a frequency of zero. In a real population, this would not occur, variants being sustained at low frequencies by a balance of mutation and purifying selection (Haldane, 1937; Haigh, 1978). To correct for this we postprocessed the initial and final frequency values from our simulations before calculating our distance, imposing a minimum minority allele frequency of 0.1%. All changes in allele frequency below this threshold were ignored, such that, for example, if a variant changed from 0.5% to 0%, this was processed after the fact so that the variant changed from 0.5% to 0.1%. The choice of threshold here is conservative; leading to a conservatively low estimate of N_{e}.
Confidence intervals
Request a detailed protocolConfidence intervals for the effective population size were calculated as the overlap of 97.5% confidence intervals for the evolutionary rates in the observed data, calculated from the regression for the real data, and estimated from the simulated statistics. The overlap of these values gives an approximate 95% confidence interval for N_{e}.
Variations in methodology
Request a detailed protocolA number of choices were made in our estimation of an effective population size. The effects of each of these choices were explored through further calculation and simulation. Results are shown in Supplementary file 1.
Approximations in the WrightFisher model
Request a detailed protocolIn the calculation to set up an initial viral population, the assignment of minority alleles to sequences becomes slow at large population sizes. Our code simulated viral genomes; a variant allele was included into the population by choosing an appropriate proportion of genomes to which the variant was assigned. For greater computational efficiency we used a pseudorandom approach for choosing genomes. Given a population size N, we generated a set P of prime numbers that were each larger than N. Given some desired allele frequency q we wish to choose qN genomes to which to assign the variant. We therefore calculated the set of numbers:
where p is a prime number sampled at random from the set P, and a is a randomly chosen primitive root of p. Given this choice of a and p, the values a^{k} (where k is an integer between one and p1) form a pseudorandom permutation of the numbers from one to p1. We constructed a set of qN genomes by choosing genomes indexed in turn by the elements of this set, beginning from k = 1, and discarding values greater than N.
To achieve calculations for population sizes larger than 10^{7} we implemented a statistical averaging method. We generated a single population of size 10^{6}, then generated 200 outcomes of a single generation of the same size, recording allele frequencies in each case. In order to simulate a value of N of size r x 10^{6} we compared the frequencies of the initial population to the mean frequencies of a random set of r outcomes. This is equivalent of simulating transmission from a population of size r x 10^{6} in which the initial population contains r copies of each of one of 10^{6} genotypes.
Phylogenetic analysis
Request a detailed protocolConsensus sequences of data were analysed using the BEAST2 software package (Bouckaert et al., 2014). Consensus sequences from each viral segment were concatenated then aligned using MUSCLE (Edgar, 2004) before performing a phylogenetic analysis on the whole genome sequence alignment. The B/Venezuela/02/2016 sequence was used to root the alignment, the haemagglutinin segment of this virus having been identified as being very close to those from the patient. Trees were generated using the HKY substitution model (Hasegawa et al., 1985). A Monte Carlo process was run for 10 million iterations, generating a consensus tree with TreeAnnotator using the first 10% of trees as burnin. Figures were made using the FigTree package (http://tree.bio.ed.ac.uk/software/figtree/).
Haplotype reconstruction
Request a detailed protocolHaplotype reconstruction was performed using multilocus polymorphism data generated by the SAMFIRE software package (Illingworth, 2016). Variant loci in the genome were identified as those at which a change in the consensus nucleotide was observed between the initial and the final consensus. The shortread data were then processed, converting reads into strings of alleles observed at these loci; a single pairedend read may describe alleles at none, one, or multiple loci. Next, these strings were combined using a combinatorial algorithm to construct a list of singlesegment haplotypes, sufficient to explain all of the observed data; no frequencies were inferred at this point. Finally, a Dirichletmultinomial model was used to infer the maximum likelihood frequencies of each haplotype given the data from each time point (Illingworth, 2015). Formally, we divided reads into sets, according to the loci at which they described alleles. A multilocus variant consists of an observation of some specific alleles at the loci in question. By way of notation, we denote by ${n}_{i}^{a}$ the number of reads in set i which describe the multilocus variant a, and denote the total number of reads in the set as N_{i}. Given a set of haplotypes with frequencies given by the elements of the vector q, we write as ${q}_{i}^{a}$ the summed frequencies of haplotypes that match each multilocus variant a in set i. For example, the haplotypes ATA and ATG would both match the multilocus variant AT describing alleles at only the first two loci. We now express a likelihood for the haplotype frequencies:
Here the parameter C describes the extent of noise in the sequence data, a lower value indicating a lower confidence in the sequence data. Haplotype reconstruction was performed by finding the maximum likelihood value of the vector of haplotype frequencies q. A value of C = 200 was chosen for the calculation, representing a conservative estimate given the prior performance of the sequencing pipeline used in this study (Illingworth et al., 2017). In contrast to previous calculations in which an evolutionary model was fitted to data (Illingworth, 2015), haplotype frequencies for each time point and for each viral segment were in this case inferred independently, with no underlying evolutionary model.
Data describing influenza A/H3N2 infection
Request a detailed protocolOur analysis of data describing longterm influenza A/H3N2 infection was performed on data from a previous study (Xue et al., 2017). As our method does not require an exceptional quality of sequencing data to calculate a rate of evolution more samples were included in our analysis than were examined in the original study. Using the codes established in the previous study, we used samples from patient W from days 0, 7, 14, 21, 28, 56, 62, 69 and 76; from patient X from days 0, 7, 14, 21, 28, 42, and 72; from patient Y from days 0, 7, 14, 21, 28, 35, 48, 56, and 70; from patient Z from days 14, 15, 20, 25, 41, 48, 55, 62, and 69. An identical procedure to that used to estimate Ne from the influenza B data was applied, calculating a rate of evolution per day from sequence data, scaling this to a rate per generation (in this case a seven hour generation time was modelled [Nobusawa and Sato, 2006]), and then running simulations to estimate N_{e}. We note that the estimates of false positive rate generated for the influenza B data were applied equally in this case, due to not having equivalent data to reestimate these values. Examining the data from patient W, our distance measurements suggested potential population structure involving the samples collected on days 62 and 69; these samples were excluded from our regression analysis.
Data availability
All sequence data is taken from previous publications, and is available from the Sequence Read Archive. Where this is sensible, raw data underlying figures has been made available in files which accompany this document.

NCBI BioProjectID PRJNA364676. Longitudinal deep sequencing of human influenza A (H3N2) from immunocompromised patients.

NCBI BioProjectID PRJNA601176. Favipiravir and zanamivir clear influenza B infection in an immunocompromised child.
References

Homologous recombination is very rare or absent in human influenza A virusJournal of Virology 82:4807–4811.https://doi.org/10.1128/JVI.0268307

BEAST 2: a software platform for bayesian evolutionary analysisPLOS Computational Biology 10:e1003537.https://doi.org/10.1371/journal.pcbi.1003537

Oseltamivirresistant novel influenza A (H1N1) virus infection in two immunosuppressed patients  Seattle, Washington, 2009MMWR. Morbidity and Mortality Weekly Report 58:893–896.

Fundamental concepts in genetics: effective population size and patterns of molecular evolution and variationNature Reviews. Genetics 10:195–205.https://doi.org/10.1038/nrg2526

MUSCLE: multiple sequence alignment with high accuracy and high throughputNucleic Acids Research 32:1792–1797.https://doi.org/10.1093/nar/gkh340

Evidence for zanamivir resistance in an immunocompromised child infected with influenza B virusThe Journal of Infectious Diseases 178:1257–1262.https://doi.org/10.1086/314440

The accumulation of deleterious genes in a populationMuller's RatchetTheoretical Population Biology 14:251–267.https://doi.org/10.1016/00405809(78)900278

A mathematical theory of natural and artificial selectionTransactions of the Cambridge Philosophical Society 23:19–41.https://doi.org/10.1017/S0305004100015176

The effect of variation of fitnessThe American Naturalist 71:337–349.https://doi.org/10.1086/280722

Intrahost emergent dynamics of oseltamivirresistant virus of pandemic influenza A (H1N1) 2009 in a fatally immunocompromised patientJournal of Infection and Chemotherapy 18:865–871.https://doi.org/10.1007/s1015601204290

Individual immune selection pressure has limited impact on seasonal influenza virus evolutionNature Ecology & Evolution 3:302–311.https://doi.org/10.1038/s415590180741x

Dating of the humanape splitting by a molecular clock of mitochondrial DNAJournal of Molecular Evolution 22:160–174.https://doi.org/10.1007/BF02101694

Identifying selection in the withinhost evolution of influenza using viral sequence dataPLOS Computational Biology 10:e1003755.https://doi.org/10.1371/journal.pcbi.1003755

Fitness inference from ShortRead data: withinhost evolution of a reassortant H5N1 influenza virusMolecular Biology and Evolution 32:3012–3026.https://doi.org/10.1093/molbev/msv171

On the effective depth of viral sequence dataVirus Evolution 3:vex030.https://doi.org/10.1093/ve/vex030

Characterization of recombinant influenza B viruses with key neuraminidase inhibitor resistance mutationsJournal of Antimicrobial Chemotherapy 55:162–169.https://doi.org/10.1093/jac/dkh528

Robust estimation of recent effective population size from number of independent origins in soft sweepsMolecular Biology and Evolution 36:2040–2052.https://doi.org/10.1093/molbev/msz081

Effective population size and population subdivision in demographically structured populationsGenetics 162:501–519.

Favipiravir and zanamivir cleared infection with influenza B in a severely immunocompromised childClinical Infectious Diseases 9:ciaa023.https://doi.org/10.1093/cid/ciaa023

Comparison of the mutation rates of human influenza A and B virusesJournal of Virology 80:3675–3678.https://doi.org/10.1128/JVI.80.7.36753678.2006

Quantifying influenza virus diversity and transmission in humansNature Genetics 48:195–200.https://doi.org/10.1038/ng.3479

Transition between stochastic evolution and deterministic evolution in the presence of selection: general theory and application to virologyMicrobiology and Molecular Biology Reviews 65:151–185.https://doi.org/10.1128/MMBR.65.1.151185.2001

Oseltamivir resistance during treatment of influenza A (H5N1) InfectionYearbook of Medicine 2006:70–71.https://doi.org/10.1016/S00843873(08)703584

Comparison of NextGeneration sequencing technologies for comprehensive assessment of FullLength hepatitis C viral genomesJournal of Clinical Microbiology 54:2470–2484.https://doi.org/10.1128/JCM.0033016

Size of population and breeding structure in relation to evolutionScience 87:430–431.
Decision letter

Aleksandra M WalczakSenior Editor; École Normale Supérieure, France

Armita NourmohammadReviewing Editor; University of Washington, United States

Georgii A BazykinReviewer; Institute for Information Transmission Problems (Kharkevich Institute), Russian Federation
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
The manuscript assesses the intrahost effective population size of influenza based on longitudinal deep sequencing data from a chronic influenza B infection. Using principles modeling and statistical approaches, the authors show that the short length of a typical influenza infection is the key limiting factor upon selection at the withinhost level. The topic is important, as it sheds light on the interplay between the two scales of selection within and betweenhost in shaping the evolution of influenza virus.
Decision letter after peer review:
Thank you for submitting your article "A large effective population size for withinhost influenza virus infection" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Georgii A. Bazykin (Reviewer #2).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
We would like to draw your attention to changes in our revision policy that we have made in response to COVID19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.
Summary:
The manuscript presents a study on withinhost population genetics of influenza virus and in particular, inference of effective population size during chronic infection in immunocompromised patients. The topic is important as it explores the interplay between the two scales of selection: withinhost and betweenhost selection that shape the evolution of influenza. Based on the analysis of sequence polymorphism, authors infer a relatively large effective population size ~10^{7} during chronic infection, in contrast to previously inferred values of ~10^{2} or less during transmission and in acute infections. All of the reviewers agree that the findings in this manuscript are interesting and a large effective population would have significant implications for efficacy of selection during withinhost evolution of influenza. However, there are still some concerns regarding methodology, interpretation and presentation of the results which we would like to see addressed.
Essential revisions:
1) Comparison between chronic and acute infections:
The authors analyzed data from chronic influenza infections and concluded that the effective population size of the virus is high, including during acute infections. For instance, the authors argue that "the observed lack of withinhost variation in typical cases of influenza can be explained by the short period of infection; the stochastic effects of genetic drift do not limit the impact of positive selection". It is however not evident that the authors' estimates of effective population size from chronic infections apply to acute infections given the exponential increase and decrease of viral load that dominate the course of acute infections. In fact, it's not clear that effective population size is even a very useful concept in this case.
Also, McCrone et al., 2018, and Xue and Bloom, 2020, have both shown that withinhost variation in acute infections is dominated by nonsynonymous mutations, and Xue and Bloom, 2020, also document stopcodon mutations within acute infections that are rarely found at appreciable frequencies in chronic infections. These observations suggest that selection is inefficient within hosts in acute infections, contrary to the authors' claims.
Moreover, McCrone et al. see radical changes in variant frequencies over the course of a few days (Figure 2E in that work) – but lineages in chronic infections (this work) persist for many months. If the authors think that N_{e} is comparable between acute and chronic infections, how do they explain the lack of diversity observed in acute infections? One way to explain this is to maintain a high N_{e} but with strong transmission bottleneck to impose stochasticity. But as point out above, "N_{e}" is really not a welldefined quantity in this case. Alternatively, could the difference imply a lower census size in acute infections, and if so, is this consistent with differences in viral load? This issue is important in view of the proposed relevance of high N_{e} for longterm influenza evolution (e.g., last phrase of the Abstract and the last phrase of the Introduction).
Overall, the authors should acknowledge the differences between acute and chronic infections, and discuss their estimates in light of the previous observations. Moreover, it may also be helpful to revise the title to indicate that the manuscript focuses on chronic infections.
2) High N_{e} is inferred from small drift and a small rate of "substitutions" (which under the authors' terminology also account for minor changes in allele frequencies). In other words, the authors are inferring a large N_{e} based on the longerterm coexistence of multiple lineages within a host. Therefore, it would be important that the manuscript also discusses alternative explanations that could lead to such patterns of polymorphism. Importantly, as N_{e} in the manuscript is inferred from a WrightFisher (WF) model, violations in the underlying assumptions of the model can bias the results. For example, one can imagine that demographic effects like population structure could be responsible for longterm coexistence and survival of lineages, e.g., if each of the samples represents a mixture of persistent subpopulations? The authors seem to suggest this by analyzing clades A and B separately, Results and Discussion, second paragraph. Alternatively, could balancing selection in the host be responsible for maintaining this polymorphism (seems unlikely, but still a formal possibility)? A discussion and/or analysis of such alternative scenarios would be useful in assessing the robustness of the manuscript's findings.
3) Robustness of the analysis and proposed statistics:
a) It would be useful to have a clearer sense of the sensitivity of N_{e} to the cutoffs used. While a lot of care has gone into the choice, some diagrams showing the sensitivity of N_{e} to cutoff choice would better demonstrate the degree to which it is a function of low frequency variants in a straightforward way.
b) To estimate how N_{e} affects changes in allele frequencies, the authors simulate a single generation of WrightFisher evolution using initial allele frequencies from a randomly selected sample from the infection. As the equation in the subsection “Summary” indicates, populations with highfrequency alleles will experience larger changes in allele frequency at a given effective population size, so the initial distribution of allele frequencies from this randomly chosen sample can have a major effect on the expected change in allele frequencies. The authors show in Figure 2—figure supplement 1 that mutations can reach frequencies of 2030% in neuraminidase, and in the influenza A patients analyzed in Figure 2—figure supplement 3, many mutations reach these and even higher frequencies, particularly at later points in the infection. The authors should run their WrightFisher simulations with different initial allele frequencies to evaluate how this choice of allele frequencies may affect estimates of effective population size.
c) The authors design statistic D to assess their estimation of N_{e}. This statistic is a sum of changes in variant frequencies across sites (subsection “Calculation of evolutionary rates”), which is then compared between data and WrightFisher simulations for different N_{e} values. The authors seem to suggest that D should be more robust to noise (subsection “Summary”), without providing any evidence. In particular, the authors should clearly state how the assumptions they made about recombination structure in WF simulation could impact the statistics D and the interpretation of the inferred N_{e}. From the manuscript it is not clear whether WF simulations are done at the sitewise, segmentwise, or genomewise level, which would impact the correlation between changes in variant frequencies. For example, simulations done with high (free) recombination would expect a lower variance D compared to the case with strong linkage (data), for the same N_{e}. These points should be better clarified.
4) In Figure 1A, it is clear (and the authors also mention) that the patient's viral load drops to undetectable levels for over a month of the infection, and viral load also varies substantially while the patient is continually infected. Effective population size and census population size are not always directly related, but the authors should discuss how changing population sizes affect their estimate of effective population size and whether a single effective population size is adequate to represent the infection.
5) The authors calculate sequence distance between every pair of sequenced timepoints to reduce the influence of noise from sequencing error, but as a result, the points in Figure 2A are nonindependent and may contribute to a tighter confidence interval around the evolutionary rate than is realistic. In particular, changes in variant frequencies that take place during the middle of the infection will be overcounted in these pairs and will disproportionately influence the overall estimate of evolutionary rates. When the authors estimate the evolutionary distance between consecutive timepoints and divide by the number of days between them, how well does the estimate correspond to the estimates in Figure 2? What is the variance in these estimates?
6) The regression performed in Figure 2A, C, and analogous figures may be especially influenced by the few points at the right end of the distribution, which represent evolutionary distances between points spaced further apart in time. How robust is the estimate of evolutionary rate to removal of these points, or by calculation of evolutionary rate as suggested in comment 4?
7) The authors chose to infer effective population size using variants and haplotypes on the neuraminidase and hemagglutinin segments. This is an odd choice since these regions tend to experience the strongest selection, which can strongly influence the estimates of effective population size. Selection can act on linked haplotypes across the genome in some cases, but have the authors tested to see if these results hold for other gene segments as well?
8) Why are the effective population size estimates for the clade B samples calculated separately from the clade A samples? It's not evident from the SAMFIRE inference of haplotypes that clades A and B constitute separate subpopulations; it seems that they could be distinct genotypes in a wellmixed population as well, as might result from a coinfection.
9) The authors assume the generation time of 10 hours per generation for influenza B. However, if generations are longer in immunocompromised individuals, the analysis would lead to an overestimation of N_{e}. Given that the main result in this manuscript is that N_{e} is high, this possibility should at least be discussed.
https://doi.org/10.7554/eLife.56915.sa1Author response
Essential revisions:
1) Comparison between chronic and acute infections:
The authors analyzed data from chronic influenza infections and concluded that the effective population size of the virus is high, including during acute infections. For instance, the authors argue that "the observed lack of withinhost variation in typical cases of influenza can be explained by the short period of infection; the stochastic effects of genetic drift do not limit the impact of positive selection". It is however not evident that the authors' estimates of effective population size from chronic infections apply to acute infections given the exponential increase and decrease of viral load that dominate the course of acute infections. In fact, it's not clear that effective population size is even a very useful concept in this case.
Also, McCrone et al., 2018, and Xue and Bloom, 2020, have both shown that withinhost variation in acute infections is dominated by nonsynonymous mutations, and Xue and Bloom, 2020, also document stopcodon mutations within acute infections that are rarely found at appreciable frequencies in chronic infections. These observations suggest that selection is inefficient within hosts in acute infections, contrary to the authors' claims.
Moreover, McCrone et al. see radical changes in variant frequencies over the course of a few days (Figure 2E in that work) – but lineages in chronic infections (this work) persist for many months. If the authors think that N_{e} is comparable between acute and chronic infections, how do they explain the lack of diversity observed in acute infections? One way to explain this is to maintain a high N_{e} but with strong transmission bottleneck to impose stochasticity. But as point out above, "N_{e}" is really not a welldefined quantity in this case. Alternatively, could the difference imply a lower census size in acute infections, and if so, is this consistent with differences in viral load? This issue is important in view of the proposed relevance of high N_{e} for longterm influenza evolution (e.g., last phrase of the Abstract and the last phrase of the Introduction).
Overall, the authors should acknowledge the differences between acute and chronic infections, and discuss their estimates in light of the previous observations. Moreover, it may also be helpful to revise the title to indicate that the manuscript focuses on chronic infections.
We acknowledge that it is important to relate our result, derived from an unusual case of infection, to more regular cases of influenza in humans. We first note that what is meant in our case by the effective population size is that statistic as it relates to an established influenza infection; our data do not describe the initial founding and growth of the viral infection.
Our primary point of reference to regular influenza infection comes via measurements of CT score relating to viral infection. Our reference on this suggests about 10 fewer units of CT, or close to a 1000fold numerical drop in census population size, in nonhospitalised, as opposed to hospitalised childhood cases. We now include the very rough calculation that this would suggest an N_{e} of around 10^{4} for such cases, cautioning that this is for an established infection, after the initial period of expansion from the transmission bottleneck.
We believe our consideration of an established population to match that of other studies of data from withinhost influenza infection; in order for data to be collected from such infections, the viral population must be of some minimal consensus size. Noting that the threshold frequency at which the effect of selection outweighs that of drift is 1/N_{e}s, we believe that within the window for which data can be collected, selection of 1% or greater per generation will dominate drift at an allele frequency of 1% or more. In this sense genetic drift does not limit positive selection.
On previous findings we do not completely recognise the statement that withinhost variation in acute infections is dominated by nonsynonymous mutations. If a simple count of variants is made, the majority will likely be nonsynonymous, however this reflects a fact that the large majority of possible variants are nonsynonymous for at least one viral protein. McCrone et al. state that their data, ‘suggest significant purifying selection within hosts’ while Xue and Bloom state that ‘synonymous mutations accumulate about twice as quickly as nonsynonymous mutations within hosts’. Synonymous mutations are relatively more common than nonsynonymous mutations at low frequencies, consistent with purifying selection.
We are not fully convinced that stop mutations are lethal in the traditional sense due to the nature of the influenza virus; the genome encapsulated within a virus (i.e. encoded in the RNA within a set of viral proteins) is not necessarily the same genome that was translated to produce the proteins. As such the observation of stop mutations at very low frequencies is not entirely inconsistent with efficient purifying selection. We are not aware of a great deal of work looking at stop mutations in chronic influenza infection, however the presence of purifying selection would again explain such a lack.
While we greatly admire the work of McCrone et al., we are not convinced that the withinhost changes in allele frequency that they observe are caused purely by genetic drift. Selection, population structure, and rare sequencing error could all contribute to the changes observed, and the data described in that paper, with two samples collected from each individual, do not allow for discrimination between drift and these other factors. In our case, where we have multiple samples from a host, we observe both large differences between individual samples (in common with McCrone) but also an underlying pattern that suggests a large withinhost population size. Previous data describing withinhost evolution does not contradict our result.
We have revised the title to, ‘A large effective population size for established influenza infection’. This recognises that our inference neglects the initial phase of viral growth, which may arise from a single particle. While we cannot with our method directly evaluate N_{e} for acute infections, we believe that an argument based on CT values carries some weight when applied to these cases.
2) High N_{e} is inferred from small drift and a small rate of "substitutions" (which under the authors' terminology also account for minor changes in allele frequencies). In other words, the authors are inferring a large N_{e} based on the longerterm coexistence of multiple lineages within a host. Therefore, it would be important that the manuscript also discusses alternative explanations that could lead to such patterns of polymorphism. Importantly, as N_{e} in the manuscript is inferred from a WrightFisher (WF) model, violations in the underlying assumptions of the model can bias the results. For example, one can imagine that demographic effects like population structure could be responsible for longterm coexistence and survival of lineages, e.g., if each of the samples represents a mixture of persistent subpopulations? The authors seem to suggest this by analyzing clades A and B separately, Results and Discussion, second paragraph. Alternatively, could balancing selection in the host be responsible for maintaining this polymorphism (seems unlikely, but still a formal possibility)? A discussion and/or analysis of such alternative scenarios would be useful in assessing the robustness of the manuscript's findings.
For additional clarity we note that our rate is equivalent to a number of substitutions per day.
Rather than the coexistence of multiple lineages we infer a rate of N_{e} based upon the rate of change within (primarily clade A) of the viral population. Multiple lineages are not required in the sense that there could be a fully wellmixed population and we could still infer N_{e} using our method. Explicitly, we derive a rate of change in the viral population, measured across multiple samples, and identify a WrightFisher population which under genetic drift matches this rate of change.
We have added a few words to clarify the cladal structure of the population. We believe that the infection is founded by a single viral population (as opposed to coinfection) and that subsequently there is a branching event, so that clades A and B become spatially separated in the host and evolve independently of one another. Our guess is that the lessfrequently observed clade B includes a smaller number of viruses, and so evolves faster under genetic drift.
We note two possible deviations from our WrightFisher model. Firstly, population structure going beyond the simple cladal structure we observe would lead to a reduction in the value of N_{e}; we cite Whitlock and Barton on this point. Such population structure would alter the value we derive i.e. it will decrease N_{e} relative to a wellmixed population, leading to an increase in the rate of change of the population that our model will detect. Regarding population structure, we note that a nonwellmixed population could lead to nonrepresentative sampling of the population and thereby increased distances between individual samples; this effect is included in our ‘error’ terminology in the method.
Secondly, we note the potential for selection to shape the population, noting the emergence of zanamivir resistance. Such selection would not be accounted for in our model i.e. to the extent that it is present it would increase the rate of change of the population that we will detect, but will attribute to a lower N_{e}. In this sense the presence of positive selection would lead us to underestimate N_{e}. Purifying selection is difficult to model; within the WrightFisher framework all selection is identical in leading to changes in allele frequencies with time. This has the consequence that as N_{e} becomes high the change in the population does not tend to zero. We note that there will be effects other than genetic drift affecting the population, and stick to our definition that our effective population size is the size at which an idealised population evolving under drift matches the behaviour of our data.
3) Robustness of the analysis and proposed statistics:
a) It would be useful to have a clearer sense of the sensitivity of N_{e} to the cutoffs used. While a lot of care has gone into the choice, some diagrams showing the sensitivity of N_{e} to cutoff choice would better demonstrate the degree to which it is a function of low frequency variants in a straightforward way.
We have made two significant cutoffs in our method. The first is to remove what we believe to be false positive variant calls in the sequence data, while the second is to impose a hard cut of 0.1% allele frequency when making our calculation of distance. To evaluate these we have rerun our calculations in a way that removes each of these in turn; we find that the resulting change in N_{e} is not greatly changed by either of these. We have added Supplementary file 1 which contains inferences for calculations run with parameters other than the default parameters.
b) To estimate how N_{e} affects changes in allele frequencies, the authors simulate a single generation of WrightFisher evolution using initial allele frequencies from a randomly selected sample from the infection. As the equation in the subsection “Summary” indicates, populations with highfrequency alleles will experience larger changes in allele frequency at a given effective population size, so the initial distribution of allele frequencies from this randomly chosen sample can have a major effect on the expected change in allele frequencies. The authors show in Figure 2—figure supplement 1 that mutations can reach frequencies of 2030% in neuraminidase, and in the influenza A patients analyzed in Figure 2—figure supplement 3, many mutations reach these and even higher frequencies, particularly at later points in the infection. The authors should run their WrightFisher simulations with different initial allele frequencies to evaluate how this choice of allele frequencies may affect estimates of effective population size.
The reviewers are correct that the allele frequencies used to initiate the WrightFisher model may affect the inferred effective population size. When we calculated replicate simulated populations we accounted for this; in each of the replicate simulations a random sample from the population was chosen to provide the allele frequencies for each segment of the simulated population. The uncertainty bars in our calculations therefore incorporate the uncertainty intrinsic to the initial choice of allele frequencies.
c) The authors design statistic D to assess their estimation of N_{e}. This statistic is a sum of changes in variant frequencies across sites (subsection “Calculation of evolutionary rates”), which is then compared between data and WrightFisher simulations for different N_{e} values. The authors seem to suggest that D should be more robust to noise (subsection “Summary”), without providing any evidence. In particular, the authors should clearly state how the assumptions they made about recombination structure in WF simulation could impact the statistics D and the interpretation of the inferred N_{e}. From the manuscript it is not clear whether WF simulations are done at the sitewise, segmentwise, or genomewise level, which would impact the correlation between changes in variant frequencies. For example, simulations done with high (free) recombination would expect a lower variance D compared to the case with strong linkage (data), for the same N_{e}. These points should be better clarified.
We have now incorporated further explanation into the Materials and methods, describing in a more formal manner how our statistic works and why it is more robust than a simple distance metric based upon pairs of samples from a population. We have clarified that our WF simulations were done at the genomewise level. Based on prior evidence from human infection, simulations assumed an absence of intrasegment recombination or of reassortment between segments; this is now more clearly stated.
4) In Figure 1A, it is clear (and the authors also mention) that the patient's viral load drops to undetectable levels for over a month of the infection, and viral load also varies substantially while the patient is continually infected. Effective population size and census population size are not always directly related, but the authors should discuss how changing population sizes affect their estimate of effective population size and whether a single effective population size is adequate to represent the infection.
The calculation we make in Clade A was performed over the samples in this clade used up to the point at which favipiravir was first used; this is shown by the green box in Figure 1A. Our belief is that CT score is somewhat noisy, sometimes providing a better measurement of the amount of viral material on a swab than of the consensus population size. Previous modelling of these data suggests a smooth, roughly 8fold decline in viral load during this period (Lumby et al., 2020). We believe that the prefavipiravir set of samples is the most appropriate one from which to derive a headline figure for effective population size, the subsequent clinical intervention being an unusual event.
We have added a note that our estimate for clade B spanned the interval in time with the bottleneck; this may be a reason for its lower value. We also note that our estimate is of a mean effective population size.
5) The authors calculate sequence distance between every pair of sequenced timepoints to reduce the influence of noise from sequencing error, but as a result, the points in Figure 2A are nonindependent and may contribute to a tighter confidence interval around the evolutionary rate than is realistic. In particular, changes in variant frequencies that take place during the middle of the infection will be overcounted in these pairs and will disproportionately influence the overall estimate of evolutionary rates. When the authors estimate the evolutionary distance between consecutive timepoints and divide by the number of days between them, how well does the estimate correspond to the estimates in Figure 2? What is the variance in these estimates?
We have provided further explanation of our method. Our basic rationale is that individual samples from the population are considerably affected by error, such that the error in each sample is larger than the true evolutionary distances undergone by the population. The raw distances between samples in consecutive timepoints are shown in the leftmost points in Figure 2A; the mean of these values is 42.6 (standard deviation 8.5), with essentially no correlation between these values and the number of days which separate the samples (pvalue 0.97 from a correlation test performed in Mathematica 11). If we persist with this calculation we obtain a mean change per day in the population of just over 9 nucleotides per day, greater than the total inferred change of close to 8 nucleotides for clade A across five months of evolution. We believe that noise in the individual samples greatly outweighs the genuine signal of evolutionary change in the population in such a way that the simple comparison of pairwise samples does not produce an accurate result.
Assuming a constant underlying effective population size (or failing that, calculating some kind of mean), our regression allows us to infer a rate of evolution even in the presence of considerable noise. We acknowledge that changes in the population during the middle of infection are overrepresented but do not have a solution to this; the use of multiple samples is intrinsic to our approach.
6) The regression performed in Figure 2A, C, and analogous figures may be especially influenced by the few points at the right end of the distribution, which represent evolutionary distances between points spaced further apart in time. How robust is the estimate of evolutionary rate to removal of these points, or by calculation of evolutionary rate as suggested in comment 4?
We have checked the slope of the regression by removing points from either the beginning and the end of the infection. In Figure 2A, removing between zero and four time points from each end of the data in any combination gives 25 regression coefficients; all of these fall within the 97.5% confidence interval that we report for our original calculation. [Note : Removing a time point removes all distances associated with that point, removing multiple points from the figure]. In Figure 2C there are data from only four time points; here removing either the first or the last leads to regression coefficients within the original confidence interval. We believe that the use of consecutive samples to assess effective population size gives a highly misleading result due to noise and confounding factors in the data greatly exaggerating the real rate of change of the population. We therefore omit this result from the main text, though we note that performing this calculation for clade A gives an effective population size of approximately 800.
7) The authors chose to infer effective population size using variants and haplotypes on the neuraminidase and hemagglutinin segments. This is an odd choice since these regions tend to experience the strongest selection, which can strongly influence the estimates of effective population size. Selection can act on linked haplotypes across the genome in some cases, but have the authors tested to see if these results hold for other gene segments as well?
This is a misunderstanding of our approach. We illustrate the cladal structure of the population using haplotype reconstructions calculated using sequence data for neuraminidase and haemagglutinin. These segments were chosen as they had slightly higher levels of genetic diversity, giving the clearest illustrations of a pattern that was visible across all of the viral segments. However, the calculation of effective population size was calculated genomewide, using data from all viral segments. We have amended the text to greater highlight the illustrative nature of the haplotype reconstructions we present.
8) Why are the effective population size estimates for the clade B samples calculated separately from the clade A samples? It's not evident from the SAMFIRE inference of haplotypes that clades A and B constitute separate subpopulations; it seems that they could be distinct genotypes in a wellmixed population as well, as might result from a coinfection.
We believe that it is unlikely that these clades arise from a wellmixed population. The samples we have are deepsequenced, generally to in excess of 2000x coverage. However, considerable differences are observed between these samples; in a wellmixed population, the samples might exhibit a pattern of evolution, but evolution would follow a continuous pattern of change rather than identifiably (in the haplotype reconstruction) being from one subpopulation or another. Our ‘clade B’ samples describe the 18^{th}, 40^{th}, and 41^{st} samples from the host. The 18^{th} sample is evolutionarily intermediate between everything in ‘clade A’ and the final two samples. This leads us to the belief that the two clades begin as a single transmitted population (i.e. not from a coinfection), that clades A and B are very largely spatially separate, and that clade B evolves away from clade A over time, potentially as a result of genetic drift. We have added further detail to the text describing the observed relationship between samples.
9) The authors assume the generation time of 10 hours per generation for influenza B. However, if generations are longer in immunocompromised individuals, the analysis would lead to an overestimation of N_{e}. Given that the main result in this manuscript is that N_{e} is high, this possibility should at least be discussed.
We are not aware of a biological reason why the generation time for influenza would be different for immunocompromised individuals, but acknowledge that this parameter might contain some uncertainty. We have explored the effect of changes in the generation time in Supplementary file 1.
https://doi.org/10.7554/eLife.56915.sa2Article and author information
Author details
Funding
Wellcome (101239/Z/13/Z)
 Christopher JR Illingworth
Wellcome (101239/Z/13/A)
 Christopher JR Illingworth
Wellcome (105365/Z/14/Z)
 Casper K Lumby
Isaac Newton Trust
 Christopher JR Illingworth
Helsingin Yliopisto
 Christopher JR Illingworth
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Senior Editor
 Aleksandra M Walczak, École Normale Supérieure, France
Reviewing Editor
 Armita Nourmohammad, University of Washington, United States
Reviewer
 Georgii A Bazykin, Institute for Information Transmission Problems (Kharkevich Institute), Russian Federation
Version history
 Received: March 13, 2020
 Accepted: July 30, 2020
 Accepted Manuscript published: August 10, 2020 (version 1)
 Version of Record published: August 17, 2020 (version 2)
Copyright
© 2020, Lumby et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,324
 Page views

 219
 Downloads

 8
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Evolutionary Biology
 Genetics and Genomics
Microbial plankton play a central role in marine biogeochemical cycles, but the timing in which abundant lineages diversified into ocean environments remains unclear. Here, we reconstructed the timeline in which major clades of bacteria and archaea colonized the ocean using a highresolution benchmarked phylogenetic tree that allows for simultaneous and direct comparison of the ages of multiple divergent lineages. Our findings show that the diversification of the most prevalent marine clades spans throughout a period of 2.2 Ga, with most clades colonizing the ocean during the last 800 million years. The oldest clades – SAR202, SAR324, Ca. Marinimicrobia, and Marine Group II – diversified around the time of the Great Oxidation Event, during which oxygen concentration increased but remained at microaerophilic levels throughout the MidProterozoic, consistent with the prevalence of some clades within these groups in oxygen minimum zones today. We found the diversification of the prevalent heterotrophic marine clades SAR11, SAR116, SAR92, SAR86, and Roseobacter as well as the Marine Group I to occur near to the Neoproterozoic Oxygenation Event (0.8–0.4 Ga). The diversification of these clades is concomitant with an overall increase of oxygen and nutrients in the ocean at this time, as well as the diversification of eukaryotic algae, consistent with the previous hypothesis that the diversification of heterotrophic bacteria is linked to the emergence of large eukaryotic phytoplankton. The youngest clades correspond to the widespread phototrophic clades Prochlorococcus, Synechococcus, and Crocosphaera, whose diversification happened after the Phanerozoic Oxidation Event (0.45–0.4 Ga), in which oxygen concentrations had already reached their modern levels in the atmosphere and the ocean. Our work clarifies the timing at which abundant lineages of bacteria and archaea colonized the ocean, thereby providing key insights into the evolutionary history of lineages that comprise the majority of prokaryotic biomass in the modern ocean.

 Evolutionary Biology
 Genetics and Genomics
In many species, meiotic recombination events tend to occur in narrow intervals of the genome, known as hotspots. In humans and mice, double strand break (DSB) hotspot locations are determined by the DNAbinding specificity of the zinc finger array of the PRDM9 protein, which is rapidly evolving at residues in contact with DNA. Previous models explained this rapid evolution in terms of the need to restore PRDM9 binding sites lost to gene conversion over time, under the assumption that more PRDM9 binding always leads to more DSBs. This assumption, however, does not align with current evidence. Recent experimental work indicates that PRDM9 binding on both homologs facilitates DSB repair, and that the absence of sufficient symmetric binding disrupts meiosis. We therefore consider an alternative hypothesis: that rapid PRDM9 evolution is driven by the need to restore symmetric binding because of its role in coupling DSB formation and efficient repair. To this end, we model the evolution of PRDM9 from first principles: from its binding dynamics to the population genetic processes that govern the evolution of the zinc finger array and its binding sites. We show that the loss of a small number of strong binding sites leads to the use of a greater number of weaker ones, resulting in a sharp reduction in symmetric binding and favoring new PRDM9 alleles that restore the use of a smaller set of strong binding sites. This decrease, in turn, drives rapid PRDM9 evolutionary turnover. Our results therefore suggest that the advantage of new PRDM9 alleles is in limiting the number of binding sites used effectively, rather than in increasing net PRDM9 binding. By extension, our model suggests that the evolutionary advantage of hotspots may have been to increase the efficiency of DSB repair and/or homolog pairing.