The sequence diversity of viral populations within individual hosts is the starting 24 material for selection and subsequent evolution of RNA viruses such as foot--and--mouth 25 disease virus (FMDV). Using next--generation sequencing (NGS) performed on a Genome 26
Analyzer platform (Illumina), this study compared the viral populations within two 27 bovine epithelial samples (foot lesions) from a single animal with the Inoculum used to 28 initiate experimental infection. Genomic sequences were determined in duplicate 29 sequencing runs, and the consensus sequence determined by NGS, for the Inoculum, 30 was identical to that previously determined using the Sanger method. However, NGS 31 reveals the fine polymorphic sub--structure of the viral population, from nucleotide 32 variants present at just below 50% frequency to those present at fractions of 1%. Some 33 of the higher frequency polymorphisms identified encoded changes within codons 34 associated with heparan sulphate binding and were present in both feet lesions 35 revealing intermediate stages in the evolution of a tissue--culture adapted virus 36
replicating within a mammalian host. We identified 2,622, 1,434 and 1,703 37 polymorphisms in the Inoculum, and in the two foot lesions respectively: most of the 38 substitutions occurred only in a small fraction of the population and represent the 39 progeny from recent cellular replication prior to onset of any selective pressures. We 40 estimated an upper limit for the genome--wide mutation rate of the virus within a cell to 41 be 7.8 x 10 --4 per nt. The greater depth of detection, achieved by NGS, demonstrates that 42 this method is a powerful and valuable tool for the dissection of FMDV populations 43 within--hosts. 44 45 Introduction 46 RNA viruses evolve rapidly due to their large population size, high replication rate and 47 poor proof--reading ability of their RNA--dependent RNA polymerase. These viruses exist 48 as heterogeneous and complex populations comprising similar but non--identical 49 genomes, but the evolutionary importance of this phenomenon remains unclear (15, 16, 50 25). Consensus sequencing identifies the predominant or major viral sequence present 51 in a sample, but is uninformative about minority variants that are present. Evidence for 52 population heterogeneity, where individual sequences differ from the consensus 53 sequence, has been routinely obtained using cloning approaches (1, 8), providing 54
insights into the evolutionary processes that shape viral populations. Unfortunately, 55 these cloning processes are laborious and usually provide only a limited resolution of 56 the mutant spectrum within a sample. 57
Next--Generation Sequencing (NGS) techniques offer an unprecedented 'step--change' 58 increase in the amount of sequence data that can be generated from a sample. Albeit 59 mostly used for de--novo sequencing of large genomes, NGS can be applied to re--60 sequence short viral genomes to obtain an ultra--deep coverage. Therefore, NGS has the 61 potential to provide information beyond the consensus for a viral sample by revealing 62 nucleotide substitutions present in only a small fraction of the population. 90  The  samples  analysed  were  collected  during  an  infection  experiment,  in  which  a  single  91  bovine  host  was  inoculated  intradermolingually  with  a  dose  of  10 5.7 50%  tissue  culture  92  infective doses (TCID50) of FMDV (O1BFS 1860). The full--length FMDV genome sequence  93  of this sample had been previously determined using Sanger sequencing (EU448369)  94 and was used as a reference genome in this study. The Inoculum was derived from a 95 bovine tongue vesicle specimen that had been passaged extensively in cell culture (9). 96
Total RNA (TRIzol, Invitrogen, Paisley, UK) was extracted from a sample of the Inoculum 97 as well as two 10% tissue suspensions prepared from epithelial lesions (front left foot 98
[FLF] and back right foot [BRF]) collected from the animal at 2 days post inoculation. 99
Reverse transcription was performed using an enzyme with high specificity 100 (Superscript III reverse transcriptase, Invitrogen), and an oligo--dT primer (see Table 1 ). 101
For each sample, two PCR reactions generating long overlapping fragments (4557 bp  102 and 4317 bp respectively) were carried out using a proof--reading enzyme mixture 103
(Platinum Taq In order to make direct comparisons between the two runs, we trimmed reads from the 138 second run to 50nt. Typically, quality scores decreased along a read, as the reliability of 139 the sequencing process decreased with the number of cycles of the Sequencing 140
Platform. The second run yielded much better qualities following an upgrade of the 141
Illumina platform. For both runs, reads with average error per nt below a fixed 142 threshold (θ = 0.2%) were discarded to generate a flatter error profile along the read 143 (see Appendix and Figure A1 ). The first and last 5 nts of each aligned read were 144 removed from the analysis as they showed a higher number of mismatches to the 145 reference sequence due to insertions or deletions close to the edges of the reads ( Figure  146 A2). More details can be found in Appendix. 147
Validation and analysis of sequence diversity in the samples 148 The frequency of site--specific polymorphisms was estimated from the frequency of 149 mismatches of the aligned reads to the reference genome. A proportion of these 150 mismatches was expected to be artifactual, arising from a base mis--calling in the 151 sequencing process, or from a PCR error in the amplification of the sample. In order to 152 identify polymorphisms arising from possible base mis--calls in the sequencing reaction, 153
we used the quality score of each nucleotide read to compute the average probability of hypothesize that the probability that PCR errors in both runs independently generated 160
identical base changes at the same site is very low. Based on values quoted for the 161 enzymes used, we estimated that the error rate for the combined RT--PCR amplification 162 process to be 7.7x10 --6 per base pair copied (2, 31, 38). We therefore defined 163 polymorphic sites that could not be attributed to sequencing errors and at which both 164 the most common and second most common nucleotides were the same between the 165 two runs to be 'qualitatively validated sites'. For each site in the set of qualitatively 166 validated polymorphisms, we computed the 95% confidence intervals for the 167 polymorphism frequency using the binomial distribution above. If the 95% confidence 168 intervals from each run overlapped, we defined the polymorphism frequency estimates 169 from the two runs to be in quantitative agreement. 170
We assessed the quantitative repeatability of site--specific polymorphism frequency 171 estimates by calculating Spearman Rank correlation coefficients between 172 polymorphism frequencies in the samples within each run and between polymorphism 173 frequencies from runs 1 and 2. 174
We counted the number of transitions (Ts) and transversions (Tv) observed at 175 qualitatively validated sites across the genome, we computed κ=2Ts/Tv, and the 176 relative distribution of mutations across the 1 st , 2 nd , and 3 rd codon positions across the 177 open reading frame (ORF). We obtained an estimate for dN/dS as follows: for each 178 codon of the reference ORF, we computed the expected number of synonymous (si) and 179 non--synonymous sites (ni) and, for each read j spanning that codon, the number of 180 observed synonymous (s r ij) and non--synonymous substitutions (n r ij). Using all the 181 codons where si>0, we then obtained an estimate of the number of synonymous 182 substitutions per synonymous site (pS) and non--synonymous substitutions per non--183 synonymous site (pN):
(and analogously for pN) ,where mi is the 184 number of reads covering codon i and ncod is the total number of codons in the ORF. 185 From pN and pS we obtained dN/dS according to (34). 186
We calculated the number of validated sites at which STOP--codons are observed within 187 the reading frame, and used these count to estimate an upper limit on the mutation rate. where N is the number of sites and piX is 218 the fraction of reads bearing nt X at site i. The entropy measures the amount of 219 "disorder" in the population, and it is maximum at a site when all four bases are equally 220
represented. 221 222
In this section, we discuss the results of the for BRF ( Figure 1A Distribution of polymorphisms across the genome 300 There were 12 SSPs, whose average frequency between the two runs is above 1% in the 301
Inoculum, 19 in FLF, and 25 in BRF (see Supplementary Table) . Some of these were 302 clustered in the capsid protein region (beginning of protein VP3) (1 in the Inoculum, 4 303 in FLF and 5 in BRF) and in the 3' untranslated region (UTR) (6 in the Inoculum, 5 in 304 FLF and 6 in BRF). Where single reads spanning these sites within the VP3 or 3'UTR 305
were available, there was no evidence that that these mutations were linked together on 306 individual FMDV genomes. In particular, the first cluster was shared between the two 307 foot samples and corresponded to changes encoding amino acid residues associated 308 with heparan sulphate (HS) binding. The Inoculum used in this experiment had 309 undergone extensive cell culture passage and, in common with other in--vitro adapted 310
hosts drives the reversion of positively charged amino acid residues at specific sites in 312 the viral capsid (20, 37). A consensus level substitution (>50%) exists within both feet 313 samples of run 1 compared to the reference sequence (see above and the  314  Supplementary Table) . This polymorphism corresponded to a change within the 60 th 315 codon of protein VP3 (VP3 60 ). Although below the level of the consensus sequence, 316
additional qualitatively validated SSPs that were present in both feet samples were 317 detected at four further sites (VP2 134 , two codon positions within VP3 56 and VP3 59 ) that 318
impact on the ability of FMDV to bind HS. All but one of the mutations clustered within 319 the 3' UTR of the three samples were located within the first four RNA--RNA pairings 320 either side of the apex of a conserved stem--loop. This structure, one of two stem--loops 321
previously predicted for FMDV and other picornaviruses (6) (33) is thought to generate 322 long--distance RNA--RNA interactions that may impact upon viral replication (40). The 323
presence of shared mutations between the two foot samples suggests a common history 324
for the viruses arising as a result of the shared route of intra--host transmission from 325 initial replication sites in the tongue to epithelial sites in the feet via the blood. 326
However an alternative explanation -that the virus is subject to a common selection 327 pressure in both sites cannot be ruled out. 328 suggesting that the absence of observed genetic variability may be due to lack of power 333 to detect it. By grouping the site--specific polymorphism frequencies into discrete bins, 334
we can examine the proportion of sites experiencing different polymorphic frequencies 335
and thereby obtain a comprehensive picture of the heterogeneity in the viral 336 populations ( Figure 5 ). Across the three samples, most sites exhibit a range of low--337
frequency SSPs between 0.01% --1%. Only a few sites showed higher frequency 338 polymorphism, and these sites were more numerous for the samples from the feet than The presence of STOP codons can be used to obtain an upper limit on the mutation rate 369
(per nucleotide per transcription event) of this virus. We hypothesise that these 370 mutations are lethal and are therefore generated in the last round of cellular replication. 371
Moreover, by assuming a replication strategy involving the minimum number of 372 copying events in the cell (the "stamping machine" strategy, see (44) of the infection across hosts (5) can be explained with this mechanism, which is made 412 more accessible to study by NGS. Moreover, the statistical characteristics of the SSPs we 413 identified (κ, dN/dS) are very similar to those found previously (8), further 414
corroborating the validity of our results. Finally, randomizations of the diversity 415 measured in the capsid region allowed us to obtain simulated clones whose 416 characteristics in terms of mutation were analogous to those found in (8). We conclude 417
that NGS data can be used to examine the nucleotide diversity of each genome position 418 at unprecedented resolution. Observing the mutant spectrum of the viral population at 419 a fine resolution will provide a more sophisticated understanding of evolutionary 420 processes shaping its variability. 421
Comparisons between the sequences recovered from the Inoculum and clinical lesions 422
provide new insights into the impact of early replication events on viral evolution 423 within a host. This study reveals that only a few sites displayed mutations present in a 424 large fraction of the population, i.e. high frequency polymorphisms (>1%), while the 425 vast majority of the polymorphisms were present at lower frequencies. We hypothesize 426 that the high frequency polymorphisms have been selected over multiple rounds of 427 replication within cells, and that the lower frequency polymorphisms most likely 428 directly reflect the high rate of mutation experienced by these viruses, as our estimate 429
of an upper limit of the genome--wide mutation rate suggests. In this study we used a cell 430 culture adapted virus (as the Inoculum) which gave us the opportunity to monitor 431 changes at specific loci associated with the HS binding site that were under selection 432 pressure during initial replication in a mammalian host. Examination of these sites 433 (collated in the Supplementary Table) identifies hot and cold spots in the HRV ORF, and some minority variants, the resolution 440
is not sufficient to observe the micro--evolutionary processes whose signature lies in 441 small fractions of the viral population (<2%). Moreover, their estimation of the 442 substitution rate during the infection is based solely on the count of the nucleotides 443 changed among those analyzed: although the value is compatible to our genome--wide, 444 mutation rate, we believe that considering the cellular process of viral replication (and 445 specifically assuming the minimum number or copying events in a cell) allows us to gain 446 a better insight of the process generating variation in the viral population and obtain a 447 more stringent upper bound. 448 Figure 5 give rise to the amount of diversity we observe. The data collected in studies like this 496
can be used for building models aimed at understanding the link between the micro--497 evolution of FMDV at the cellular scale with the population heterogeneity at the host 498 scale. We anticipate that a model of viral replication across several cell generations 499 within a host will produce a more stringent upper bound to the genome--wide mutation 500 rate. 501
Although further work is required, these findings strongly suggest that data generated 502 through the use of this methodology can provide novel insights into viral evolutionary 503 dynamics at a greater resolution than previously achieved for a positive--stranded virus 504 such as FMDV. In particular, the genome wide assessment of polymorphic frequencies is 505 likely to be an important asset in the parameterization of models that can evaluate the 506 role of quasi--species dynamics in RNA virus evolution. 507 508 The quality scores associated with each nucleotide were lower on the first run and 539 decreased towards the end of reads ( Figure  A1 ). In order to make direct comparisons 540 between the two runs, we trimmed reads from the second run to 50nt. Typically, quality 541 scores decreased along a read, as the reliability of the sequencing process decreased 542 with the number of cycles of the Sequencing Platform. As Figure A1 shows, a trade off is 543 present between the number of reads kept and their quality. For both runs, we 544 discarded reads with average error per nt below θ = 0.2%, resulting in a flatter error 545 profile along the read. 546
With this choice of the threshold, 66% of the reads were retained from the first run (a  547  total Reads alignment and trimming 553 The vast majority of the filtered 50nt--reads aligned unambiguously with less than 5 554 mismatches to the reference inoculum genome, previously established using 555
conventional Sanger sequencing (9) (Genbank accession no. EU448369), (run 1: 92.5% 556
for Inoculum, 98.9% for FLF, and 97.8% for BRF, run 2: 95.8% for Inoculum, 98.4% for 557 FLF and 96.2% for BRF). The remaining reads were either ambiguously aligned reads or 558 contained a large number (>4) of mismatches to the reference sequence, and were 559 discarded from the analysis. For each sample, an almost equal number of reads were 560 derived from positive and negative strands of the viral cDNA. 561
Further filtering of the data was undertaken after alignment of the reads. Within the 562 aligned reads, mismatches occurred with similar frequency at each of the 50nt of the 563 reads, except from the edges, where a higher number of mismatches was observed 564 ( Figure A2 ). Presumably, these peaks were due to a small number of sequences with 565 insertions or deletions close to the ends of the reads: for subsequent analysis we 566 trimmed away the first and last 5 nts of each aligned read, leaving only the 40 central 567 nucleotides where the mismatch curve was flat. 568 569 