Deep Sequencing of Viral Genomes Provides Insight into the Evolution and Pathogenesis of Varicella Zoster Virus and Its Vaccine in Humans

Optimization of Targeted Enrichment

Archived diagnostic samples of vesicular fluid contain small amounts (20–150 ng) of extracted total DNA and required whole-genome amplification (WGA) using the Φ29 polymerase to generate sufficient starting material (2–3 µg) for the SureSelect (SS) custom target capture methodology (Depledge et al. 2011). Illumina sequencing was performed on 20 vaccine rash samples collected between 1988 and 2007 and three distinct batches of the VariVax (Merck) VZV vaccine from 2008, 2010, and 2012 (table 1). To determine whether any bias was associated with these methods, we prepared five sequence libraries (WGA only, SS only, an SS only replicate, WGA + SS, and direct sequencing without WGA or SS; supplementary table S1, Supplementary Material online) from a single batch of VZV vaccine. The consensus sequences for each sample were identical, while the correlation of allele frequencies between libraries was high (Pearson’s product moment R² > 0.93) at sites where total read depth was over 50 reads per base and where the variant allele was present on two or more independent reads mapping to opposite strands (fig. 1). Deviations in the variant allele frequency of up to 3.2% were observed but these were limited to comparisons between enriched and nonenriched libraries (fig. 1). This is explained by the fact that read depths in nonenriched libraries are approximately 20-fold. A comparison of replicate samples that were independently enriched showed the standard deviation from the mean allele frequency, that is, the error of the method to be 0.654%.

Viral DNA was enriched from all vaccine and vaccine-rash samples with a median of 87.3% on-target reads (table 1). Read data aligned against the pOKA reference genome (Gomi et al. 2002) showed a bias against enrichment in GC-rich regions of the variable region R1–R5 and short region genomic repeated elements (supplementary fig. S1, Supplementary Material online). Mean read depths were consistently higher than 1,000 reads per base for all samples except the nonenriched vaccine. Genome coverage exceeded 97%, except in two cases; samples O27 (90.5%) and K48 (83.8%). The paired-end sequence reads were too short to adequately map the variable R1–R5 regions and for the majority of the rashes there was insufficient material to fill these gaps using PCR and Sanger methods. Sequencing of vaccine rash samples on the Illumina GAIIx platform gave significantly more low frequency polymorphic sites (defined as having a variant allele at <10%) in vaccine rash samples when compared with samples sequenced on an Illumina MiSeq (supplementary fig. S2, Supplementary Material online). Resequencing of three samples for which sufficient DNA was available (K11, L53, and N13) on both platforms showed the polymorphisms present in the MiSeq data set to be a subset of those present in the GAIIx data set with good correlation between like data sets (K11 − r² = 0.897, N13 − r² = 0.908, and L53 − r² = 0.897; supplementary fig. S3, Supplementary Material online). In all cases, the excess of polymorphic sites in data sets from the GAIIx did not alter the consensus sequences for these samples but do confirm that a higher sequencing error rate is associated with the GAIIx, a finding that has been previously reported (Minoche et al. 2011; Loman et al. 2012; Quail et al. 2012).

Consensus Sequence Analyses

Animal alphaherpesvirus vaccines have been shown to recombine with circulating wild type and other attenuated vaccine viruses with the potential for restoration of viral virulence (Lee et al. 2012). To exclude this in vOka-based viruses recovered from rashes, Bayesian (fig. 2) and Maximum-likelihood trees (not shown) of vaccine and rash consensus sequences were constructed (Barrett-Muir et al. 2003; Loparev et al. 2007; Sengupta et al. 2007). We found no significant evidence for the vaccine rashes representing recombination with wild-type strains. BaTs analysis of the consensus phylogeny by year of vaccination excluded clustering of rash sequences by, mode of sample preparation (AI P = 0.36, MS P = 1.0, MC [high passage in vitro culture] P = 1.0, MC [low passage culture] P = 0.45, and MC [no culture] P = 1.0) or disease presentation (zoster or varicella: AI P = 0.85, PS P = 0.99, MC [varicella] P = 1.0, and MC [zoster] P = 0.64). The last indicates that viruses that establish latency and reactivate are not limited to a subset of haplotypes within the population of vOka variants that cause rashes.

Conservation of Polymorphic Sites in Vaccine

All three sequenced vaccine batches were heterogeneous with between 235 and 336 polymorphic sites being identified. Of these, 224 were present in two or more vaccine preparations (supplementary table S2, Supplementary Material online) and 207 were shared between all three vaccine preparations, of which 36 had been previously identified by Sanger or pyrosequencing (Gomi et al. 2002; Loparev et al. 2007; Quinlivan et al. 2007; Tillieux et al. 2008; Kanda et al. 2011). These 25–120 single nucleotide polymorphisms (SNPs) were limited to single vaccine batches of which 90% were low frequency (<10%). Five sites (12694 and 12779 [open reading frame {ORF} 10], 31732 [ORF 21], 82225 [ORF 42], and 106710 [ORF 62]) previously reported as polymorphic in Merck vaccine batches were fixed for the wild-type allele in all batches sequenced here (Gomi et al. 2002; Loparev et al. 2007; Quinlivan et al. 2007; Tillieux et al. 2008). The distribution of nonsynonymous, synonymous, and noncoding mutations is shown in figure 3. The vaccine allele frequencies at each of the 224 sites shared between vaccines were highly conserved (Pearson’s product moment R² > 0.98), indicating minimal batch to batch variation at vaccine loci despite each having been issued in different years and in one case a different country (table 1; supplementary fig. S4A–C, Supplementary Material online). Therefore, we infer that the vaccines sequenced in this study are genetically compared with the batches used to inoculate the subjects whose rash viruses we analyze here. This proposal is also supported by the high correlation (Pearson’s product moment R² = 0.83) between the mean allele frequencies in the vaccines and rash samples (supplementary fig. S4D, Supplementary Material online).

Vaccine Rashes Are Mono- or Oligomorphic

Because of the increased proportion of erroneous polymorphic sites introduced by GAIIx sequencing, only the variant data from rashes sequenced on the MiSeq were analyzed further (i.e., GAIIx data were only included in consensus level analyses). Between 32 and 112 polymorphisms per virus were identified in these six uncultured vaccine rashes, when compared with 235–336 in the vaccine batches. The lower genetic diversity of rashes can also be quantified by the proportion of polymorphic sites, being 0.19% for the vaccine, which is 2–4 times greater than rash diversity (0.03–0.09%). In vitro tissue culture further reduced rash virus diversity even at ≤3 passages (supplementary fig. S5, Supplementary Material online). Of the polymorphisms in the vaccine rashes, 37–93% (mean 55%) was at the sites of 224 vaccine SNPs (loci polymorphic in two or more vaccine samples). Vaccine-allele frequencies at the majority of these sites were lower than in the original vOka vaccines (supplementary table S2 and fig. S6, Supplementary Material online). These data are consistent with the conclusion that there exist one or more bottlenecks after inoculation of vaccine and the formation of rashes. There is a low level of polymorphism at novel loci which can be attributed to mutation because of the bottleneck. The rashes sequenced here were fixed for the vaccine allele at four amino acid coding positions in ORF 62 (105705, 106262, 107252, and 108111) (supplementary table S2, Supplementary Material online) and their presence may explain the mild presentation of most vaccine rashes compared with wild-type rashes (Tsolia et al. 1990; Cohen 2007). In addition to the 224 vaccine SNPs, we identified a further 39 polymorphic sites present in at least one vaccine and one vaccine rash and two present in two or more vaccine rashes but absent in the vaccines(supplementary table S3, Supplementary Material online). As we know that none of the rashes arose from the same VariVax vaccine batch (table 1), the most likely explanation is that for these 41 sites, the vaccine allele is present in the original vaccine at frequencies that are at the limit of detection by deep sequencing and therefore are not always detected as polymorphic.

De Novo Mutations Occurring after Inoculation

A total of 33 de novo mutations were identified in the consensus sequences of all rashes. Sixteen of these were nonsynonymous and occurred in eight rashes (VR1, VR2, A182B, L53, O27, T25, T17, and T61; table 2). The fixed mutations in rash A182B, a highly passaged isolate, for which plenty of DNA was available, were confirmed by Sanger sequencing. In the six rashes sequenced by the MiSeq method (supplementary table S4, Supplementary Material online), 67% of the new mutations that were not fixed were largely present at allele frequencies of 10% or less and the resequencing of K11, N13, and L53 confirmed of 62/68 (91.2%) of the new mutations in these rashes.

Selection for Rash-Forming Alleles

We have shown the segregating allele frequencies in the Merck vOka vaccine are sufficiently nonvariable between batches, that most of the vaccine alleles can be considered to have been present in the particular vaccine batch inoculated into patients (supplementary fig. S5, Supplementary Material online). We carried out analyses to evaluate whether there was any consistent change in allele frequencies between vaccine and rashes. This analysis had to take into account the dramatic loss of genetic diversity in each rash, because it involves large allele-frequency changes from the vaccine. These large differences cause the average frequencies across all rashes to deviate from the vaccine frequencies: A pattern that could be misinterpreted as evidence of consistent trends, unless appropriate compensation is made. We used two approaches to allow for this effect: A generalized linear mixed model, and a generalized linear model with quasi-binomial error. Both approaches detected a trend for the wild-type allele to increase in frequency in the rashes (P < 0.001), the mixed effect model suggesting a significant effect only at the coding sites (P < 0.001 for nonsynonymous and synonymous sites), indicative of selection.

A second analysis was used to identify SNPs with the clearest evidence for selection: those at which the vaccine allele repeatedly changed in frequency, which could indicate loci with a role in rash formation. Of the 224 vaccine polymorphic sites, there was statistical power (≥80%) to perform a binomial analysis for 7: 59591 (ORF31), 94167 (ORF 54), 102002 (intergenic), 105356 (ORF 62), 106001 (ORF 62), 107797 (ORF 62), and 108838 (ORF 62). At six of these sites, five nonsynonymous and one noncoding, the vaccine allele decreased significantly in the rashes. The vaccine allele increased significantly at the remaining synonymous site (94167 [ORF 54]) (supplementary table S5, Supplementary Material online). Using rash genotyping data from the literature and results from SNP typing of material from 36 additional rashes, we confirmed there was selection in favor of the ancestral allele for three of the five nonsynonymous sites, 105356 and 107797 (both ORF 62) which have previously been reported (Quinlivan et al. 2007), and 106001 (ORF 62; fig. 4). We analyzed a further nine nonsynonymous, three noncoding and one synonymous sites (560 [ORF 0], 19036 [ORF 13], 58595 [ORF 31], 87306 [ORF 50], 89734 [ORF 51], 90318 [ORF 51], 97479 [ORF 55], 97748 [ORF 55], 105169 [intergenic], 107599 [ORF 62], 109137 [intergenic], 109200 [intergenic], 115269 [intergenic]) for which SNP typing data from the literature and/or genotyping of a further 36 rashes provided enough power to confirm a significant change in allele frequencies (supplementary table S5, Supplementary Material online). Of these, there was a significant decrease in the frequency of the vaccine allele at six nonsynonymous sites: 560 (ORF 0), 19036 (ORF 13), 58595 (ORF 31), 97748 (ORF 55), 97479 (ORF 55), and 105599 (ORF 62); two noncoding: 109137 (ORF 62 promoter) and 105169 (61/62 intergenic); and a significant increase in the vaccine allele at one synonymous site: 89734 in ORF 51 (fig. 4; supplementary table S5 [Supplementary Material online]). In addition to positions 105169, 105356, and 107797 (all ORF 62), which have previously been reported (Quinlivan et al. 2007), three of the loci under selection, sites 106001 (K1045E) in ORF 62 domain IV, 19036 (G199R) in ORF 13, and 97479 (V495A) in ORF 55 have not previously been identified as polymorphic in the vaccine. None of the nonsynonymous loci were associated with changes in silico in known or predicted HLA class 1 epitopes (data not shown).

To evaluate the severity of the population bottleneck, we provisionally attributed alleles with a frequency of less than 1% in rashes as postbottleneck mutations. Having discarded these allele frequencies, we modeled the size of the putative population bottleneck occurring in the vaccine rash after vaccination by fitting a beta-binomial distribution around the expectation of the vaccine (estimated from the mean allele frequencies of the three vaccines). One of the fitted parameters of the beta-binomial analysis represents the probability F that two alleles randomly drawn from within a population are identical by descent. The estimated F was 0.79 (0.77–0.81, 95% CI). Matching the allele frequency spectrum, and this value of F to the output from a stochastic model of ancestry under exponential growth (we found the best fitting growth rate corresponding to the observed F, using the MLE function of R (growth rate r = 0.0703, over 202 generations), we could explain the data with a model (supplementary fig. S7, Supplementary Material online) in which most rashes were likely to be descended from one or two genotypes (P = 0.41 and 0.51, respectively) and in which the majority of mutations (97%) occurred after bottleneck.

The availability of material from vOka varicella and zoster rashes together with the VariVax vaccine strain ancestral to both provided a unique opportunity to determine whether vOka rashes caused by viruses reactivating from latency differed genetically from those causing skin rashes directly, that is, within 42 days following inoculation (fig. 5). We found no significant differences in allele frequencies between viruses directly causing rash after vaccination (fig. 5, routes A and B) and those that established latency before reactivating to cause a zoster-rash (fig. 5, route C). The data presented here reflect only studies of the VariVax vaccine produced by Merck. Although all licensed vaccines are derived from the same attenuated vOka strain, the Biken and GSK vaccines have undergone different manufacturing processes. For example, the GSK vaccine was plaque purified and is known to have a higher percentage of vaccine SNPs overall (Kanda et al. 2011). Although it is reasonable to assume that the findings reported here apply to vaccine rashes derived from other manufacturer’s vaccines, further sequencing is required to confirm this.

Repository of Data

Consensus sequences for all vaccines and vaccine rash samples sequenced in this study have been uploaded to the NR database on Genbank under accession numbers KF558371-91.

Sample Collection and Ethics

Diagnostic samples from patients with confirmed vaccine-associated rashes (both varicella- [n = 8] and zoster-like [n = 12]) collected in the UK and United States between 1988 and 2010, were retrieved from the Breuer lab cryobank. The clinical specimens, collected as part of standard clinical practice to genotype vaccine rashes occurring in recipients of vOka vaccine, were independently obtained from patients who developed skin rashes following immunization with the Merck vOka vaccine. All rashes were confirmed as vOka strain using previously published methods (Larussa et al. 1998; Parker et al. 2006). Rashes occurring within 42 days of vaccination were classified as varicella like. The exception was one case, T61, which was reported as a varicella-like rash 90 days after vaccination. Although by our definition this was thought to be more likely to be disseminated zoster, the case was included as varicella like. Cases presenting with the classical dermatomal rash were diagnosed clinically as herpes zoster. Samples were stripped of all personal identifiers other than details of vaccination and the type of rash, no patient information was available. The use of these specimens for research was approved by the East London and City Health Authority Research Ethics Committee (P/96/046: Molecular typing of cases of VZV).

In addition, three batches of the VariVax live-attenuated VZV Vaccine preparation were collected in 2008 (VVAG [USA]), 2010 (VV10 [UK]), and 2012 (VV12 [UK]) and batch and lot numbers used to verify them as independent preparations (table 1).

Cell Culture

Cultured isolates (table 1) were passaged in melanoma cells (MeWo) that were propagated in culture media (MEM [Sigma] supplemented with 10% FBS and 1% nonessential amino acids). Viruses were harvested for DNA extraction 3–4 days after inoculation when syncytia was visible.

Library Construction, Targeted Enrichment, and Sequencing

Total DNA was extracted from each sample using the QiaAMP DNA mini kit (QIAGEN) according to manufacturer’s instructions. DNA quantification was performed using a NanoDrop spectrophotometer and those with 260/280 ratios outside the range 1.7–2.1 and 260/230 ratios out the range 1.8–2.2 were further purified using the Zymoclean Genomic DNA Clean & Concentrator (Zymo Research Corp.). WGA using GenomiPhi V2 (GE Healthcare) was performed using 10 ng of starting material where <50 ng total DNA was available (table 1). Libraries were constructed as per the standard SS XT v1.3 (Agilent) and NEBNext (New England Biolabs) protocols, the latter modified to include a second PCR amplification step (6 cycles) to enable multiplexing using standard Illumina barcodes. Enrichment for VZV sequences was performed as described previously (Depledge et al. 2011). Samples were sequenced across several Illumina platforms (GAIIx, HiSeq, and MiSeq) according to availability. Specific sample preparation and sequencing metrics are shown in table 1. Concerns about bias introduced by SS and WGA methodologies were assuaged by resequencing the vaccine VV12 under multiple conditions (supplementary table S2, Supplementary Material online) as in all cases, the R² correlation between allele frequencies exceeded 0.93, and no identifiable bias was observed regardless of sample treatment.

Genome Assembly and Variant Calling

Each data set was parsed through QUASR (Watson et al. 2013) for duplicate removal and read-trimming (-q 30, -l 50) and subsequently aligned against pOka (AB097933.1) using BWA (Li and Durbin 2009). Resulting alignments were processed using SAMTools (Li et al. 2009) to generate pileup files for each sample. A consensus sequence for each data set was called with the QUASR module pileupConsensus and a 50% frequency threshold (i.e., no ambiguities were included). Variant profiling for each data set was performed using VarScan v2.2.11 (Koboldt et al. 2012) with the following parameters: basecall quality ≥20, read depth ≥50, independent reads supporting minor allele ≥2 per strand. In addition, variant calls showing directional strand bias ≥0.85 were excluded from further analyses. Consensus sequences were generated for each rash sample but iterative repeat regions (ORIS R1, R2, R3, R4, and R5) and the terminal repeat region were trimmed before analysis.

Consensus Sequence Analyses

DNA sequences were aligned using the program Mafft, v6 (Katoh 2002) with alignments checked manually; no insertions or deletions were inferred from the alignment. Prior to further analysis the VZV repeat regions R1, R2, R3, R4, R5, and OriS (Davison and Scott 1986; Hondo and Yogo 1988) were removed as these could not be accurately mapped given the limited DNA available. Bayesian phylogenetic trees were inferred using the program Beast v1.7.2 (Drummond et al. 2012). Two independent Monte Carlo-Markov chains were run for 50,000,000 iterations for the whole genome (minus repeat regions) with a thinning of 1,000. The most appropriate model of nucleotide substitution was selected using the program jModelTest (Posada 2008), which identified a general-time-reversible model of nucleotide substitution with a gamma-shaped rate distribution and a proportion of invariable sites (GTR + γ + I) as being the best site model. The phylogeny was reconstructed under a number of different tree priors including the Bayesian skyline, constant, and exponential growth models. Runs were checked for convergence and that an effective sample size of at least 200 had been achieved for all parameters. Runs were combined using LogCombiner v1.6.1, and TreeAnnotator v1.6.1 was used to obtain the highest clade credibility tree and posterior probabilities per node. Maximum likelihood phylogenies were also inferred using Mega v5.05 (Tamura et al. 2011) under the same model of evolution. All trees were visualized and drawn using Figtree v1.3.1 (http://tree.bio.ed.ac.uk/software/figtree/, last accessed September 17, 2013). We tested for evidence of recombination using the programs GARD and SBP (Pond et al. 2006), and the BootScan method contained in the RDP3 program (Martin et al. 2010). The program BaTS (Parker et al. 2008) was used to test for an association between the Bayes-inferred phylogeny and the vaccine rash phenotype (varicella-like or zoster-like). BaTS estimates the association using the association index (Wang et al. 2001), parsimony score (Slatkin and Maddison 1989), and maximum monophyletic clade statistics.

Statistical Analyses

To compensate for intrinsic differences in the sequencing technologies used, their respective error profiles and fold differences between the mean read depths of individual samples, we identified a minimal criterion for defining biallelic sites that segregate our samples 1) total read depth ≥ 50, 2) the vaccine allele must appear ≥2 independent reads mapping to each strand, 3) the vaccine frequency must be ≥1%, and 4) the segregating site must appear in ≥2 vaccines and/or vaccine rash samples at the earlier-mentioned criteria. These criteria were applied only to data obtained using the MiSeq and HiSeq platforms and yielded a total of 368 sites that were segregating across the vaccines and vaccine rashes. The allele frequencies at these sites in the GAIIx-derived data were subsequently included if the vaccine allele was a match (i.e., not a sequencing error). The frequency of the vaccine allele in vaccine rash samples at the majority of these sites was effectively binary, either 0% or 100%. 224 sites are polymorphic in two or more vaccine batches while sites that were polymorphic in just the one vaccine were excluded from further analysis. Subsequently, we calculated the mean vaccine allele frequency for each of these polymorphic sites in the three sequenced vaccines, VV10, VV12, and VVAG. We investigated any consistent changes in allele frequency from the vaccine to the vaccine rashes. The variance in vaccine rash allele frequencies appeared higher than expected and we modeled this over dispersion around the expectation (the rash frequency) using two approaches: a generalized linear mixed model and a generalized linear model with quasi-binomial error. With the rash allele frequencies as the response variable, we compared four models in the generalized linear model: no intercept, intercept only, coding/noncoding sites, and synonymous/nonsynonymous/intergenic region.

The majority of the 224 vaccine polymorphic sites in the vaccine rashes were effectively binary, that is, fixed for either the vaccine or the wild-type allele. Consequently, to identify any specific sites where the change in allele frequency cannot be explained by genetic drift alone, we compared the proportion of vaccine rashes with the vaccine allele to the expected proportion (vaccine mean). We estimated the effect size (based upon the mean vaccine allele frequency in the vaccines and the proportion of vaccine rashes fixed for the vaccine allele) and power (Cohen 1988), and identified those sites for which there was ≥80% power to compare the vaccine rashes with the vaccine. Finally, we applied two-tailed binomial exact tests to those sites with ≥80% power with the probability of “success” equal to the mean vaccine allele frequency in the vaccines. We corrected for multiple analyses using FDR (Hochberg and Benjamini 1990).

We also used the power analyses to identify a set of the vaccine polymorphic sites that could be analyzed with SNP typing of an additional 36 vaccine rashes, which were readily available and/or the inclusion of previously published data (Loparev et al. 2007; Quinlivan et al. 2007; Thiele et al. 2011). Using binomial exact tests, we compared the proportion of the vaccine allele at these sites in all of the vaccine rashes (sequenced, SNP typed, and previously published) with the mean frequency in the sequenced vaccines (with a FDR correction for multiple analyses). In addition, we reanalyzed those sites that we showed to have sufficient power from the sequencing (discussed earlier) with the addition of SNP typed and previously published data. Finally, we used generalized linear modeling with binomial error to identify any differences between varicella and zoster vaccine rashes at the 224 vaccine polymorphic sites. All binomial tests and the generalized linear modeling were performed in R (www.r-project.org/index.html, last accessed September 17, 2013). Effect size and power calculations were done using the “pwr” package in R, and the modeling of rash allele frequencies used the lme4 and VGAM packages.

To infer the size of the putative population bottleneck occurring in the vaccine rash after vaccination, we fitted a beta-binomial distribution around the expectation of the vaccine (estimated from the mean allele frequencies of the three vaccines). It has been shown that the distribution of allele frequencies (assuming that loci are biallelic) under genetic drift can be modeled as a beta-binomial distribution (Balding 2003). The probability that two alleles randomly drawn are identical by descent, F, is directly related to the correlation parameter of the beta distribution. We fitted the beta-binomial using the R package VGAM. Finally, we modeled the rash sample as the product of a bottleneck followed by exponential growth. We found the best fitting growth rate corresponding to the observed F, using the MLE function of R. The expected value of F was calculated in a model of postbottleneck exponential growth to a nominal end point of million-fold increase. F was calculated as 1 − Π(1-1/N_i), where N_i is the population size in generation i. We then used a Monte Carlo simulation of ancestry with these population sizes, by artificially drawing the ancestors of a sample back through the 201 preceding generations. The sample size matched our average coverage (1,942 reads). From 1,000 such possible ancestries, we generated the allele frequency spectra that would occur for mutations occurring post bottleneck and prebottleneck. The mutations were uniformly distributed across all ancestors (in each category). We matched the allele frequency spectrum of alleles that were not detectable in the vaccines, presumably new mutations occurring postbottleneck, and compared F_IS with the output from a stochastic model of ancestry with exponential growth (implemented in R).

Additional SNP typing

SNP typing at positions 59591, 90318, 94167, 97479, 97748, and 106001 was conducted on an additional 36 vaccine rash samples derived from a range of patients across the UK and the United States. PCRs comprising 30 amplification cycles were performed in 12.5 µl reactions using either the Extensor 2x Master Mix (Thermo Fisher Scientific, Hanover Park, IL) or the Takara Ex Taq Kit (Clontech Laboratories, Mountain View, CA). Each reaction of 2 µl was run on 1% agarose gel and visualized on a Biorad Molecular Imager Gel Doc XR +Imagine System (Bio-Rad Laboratories, CA). PCR products were cleaned up using ExoSap-It from USB (Affymetric, Cleveland, OH) according to manufacturer’s specifications and sequenced on the 3130X/Genetic Analyzer from ABI (Life technologies, Foster City, CA). The sequence traces were analyzed with the Sequencher program from Gene Codes Corporation (Ann Arbor, MI). All primers used for SNP typing are shown in supplementary table S6 (Supplementary Material online).

Selection of Low-Level De Novo Variants

Mutations acquired postinoculation (i.e., de novo) were identified in data sets derived from the MiSeq under the following criteria; the variant allele occurred in ≥2 reads on opposite strands with a read depth ≥ 50 and a frequency exceeding 1%. To identify whether there was any selection operating upon the new mutations in the vaccine rashes, we used the Kolmogorov–Smirnov test to compare the total number of nonsynonymous mutations with the sum of the synonymous and noncoding mutations in all of the rashes.

In Silico Epitope Prediction

As fitness selection is often linked to immune escape, the nonsynonymous sites under selection 59591, 97479, 97748, 105356, 106001, 107599, and 107797 were examined using Metaserver (MacNamara and Kadolsky 2009) and Epipred (Heckerman et al. 2006) to identify HLA class 1 epitopes mapping to the vaccine and wild-type peptides and to calculate the predicted binding scores. No significant hits were reported for any site.