Introduction Hepatitis C virus (HCV) infection is a major cause of chronic liver disease, resulting in substantial morbidity and mortality worldwide, with between 123 and 170 million persons infected [1]. Transmission is predominantly via blood-to-blood contact associated with contaminated injection devices. Infection persists in approximately 70% of acute cases, leading to chronic hepatitis, and ultimately cirrhosis and the associated complications of liver failure and hepatocellular carcinoma [2]. The outcome of primary HCV infection is driven by the interplay between rapid viral evolution and host adaptive immune responses [3], [4]. Analogous to other RNA virus infections, development of effective vaccines and antiviral treatments has been constrained by the ability of these viruses to swiftly overcome evolutionary pressures such as host immunity and antiviral drugs [5]–[7]. The rapid rate of evolution in RNA viruses is driven by a highly error-prone RNA-dependent RNA polymerase (RdRp). For HCV, the estimated mutation rate is 1.2×10−4 substitutions per site per infected cell [8], and with a half-life of 3–5 hrs, it is estimated that at least 1012 particles are generated per day [9]. With these rapid kinetics, at least 109 variants with single- and double-nucleotide changes are likely to arise in each individual multiple times daily [10]. However, the observed viral complexity is much less than this prediction, largely due to reduced fitness of mutated variants [11]. While the error-prone viral replicase results in a continuous supply of new viral variants (genetic drift), purifying selection generated from host immune pressure and viral fitness results in culling and preferential selection of certain variants. When selection overbalances drift, a genetic bottleneck occurs (i.e. evolutionary events resulting in a reduction in genetic variation due to extinction of a significant proportion of the viral variants). Very little quantification of the within-host evolutionary dynamics of human RNA viruses has been reported, and limited information exists from experimental evolution systems [12], and from viruses infecting animals. For RNA viruses, the effective population size - broadly equated with the number of variants that will contribute genes to the next generation, a key parameter in viral evolution - is smaller than the census size. This discrepancy is often seen at the epidemiological level as result of a strong transmission bottleneck. Well-characterized within-host genetic bottlenecks have been observed upon transmission of a number of RNA viruses, including HIV [13]–[15]. A low effective population size has important consequences for viral evolution, as at small effective population numbers (Ne ), random processes (i.e. genetic drift) predominate over deterministic ones (i.e. selection) [16]. Genetic bottlenecks severely limit viral diversity and potentially limit replicative fitness of the resultant virus (reviewed in [12]). After a bottleneck event, RNA viruses undergo rapid evolution, leading to accumulation of deleterious mutations and the occurrence of rare fit variants, which can rapidly reach fixation and dominate the next population. Following a genetic bottleneck, within-host evolution of rapidly mutating viruses can be characterized by selective sweeps (i.e. the reduction or elimination of variants in the viral population as the result of strong selection pressures), with only a few variants emerging to dominate the future population. This phenomenon is common at the host population level, such as in influenza [17], and has been documented within-host for HCV [18], [19] and HIV [20]. The extent to which positive selection and random genetic drift contribute to the within-host evolution of HCV, and whether selective sweeps observed late in primary infection are related to genetic bottleneck events, remain to be resolved. A major challenge in studies of the within-host evolution of RNA viruses lies in the capacity to detect low frequency viral variants. Standard cloning techniques, and more recently single genome amplification, followed by sequencing have been used to detect variants present at frequencies as low as 10–15% [15], [21], [22]. Next generation sequencing (NGS), despite high technical errors [23], allows detection of rare variants present at less than 1% of the population [24]–[33]. Data analysis tools have now improved the ability to differentiate true biological variation from technical error [34], [35], and also enabled reconstruction of genomic regions of individual variants from short NGS reads [34], [36]. The early evolution of HCV has been investigated to date only in limited case series [18], [30], [35], [37]–[40], and with a strong bias towards symptomatic cases, who represent a small minority of those with acute infection, and who are known to have an increased likelihood of clearance [3], [4]. Few studies have investigated early infection in a longitudinal fashion [30], [35], [39], [41]. In these cases the frequency, distribution, and timing of viral mutations has been characterized only in limited regions of the virus and with low sensitivity in the detection of rare variants. Here we report a comprehensive analysis of longitudinal collected samples from four subjects identified very early in asymptomatic acute HCV infection. The aim was to quantify the number of successfully transmitted founder viruses, and to characterize the diversity and complexity of viral population across the genome over the course of the primary infection leading to clearance or chronicity. Results Subjects Four newly viremic seronegative subjects were studied (Figure 1, Table 1, Table S1). Subjects were enrolled in the Hepatitis C Incidence and Transmission Study (HITS) cohort, and had tested negative for HCV antibodies and RNA within 3–6 months prior [42], [43]. The estimated days post-infection (DPI) at enrolment ranged from 30–45 days (Table 1). Two subjects cleared the infection (686_Cl, 360_Cl), and the other two progressed to chronic infection (23_Ch, 240_Ch). A single viral genotype (GT) was detected in each subject: two with GT1a (686_Cl, and 23_Ch), and two with GT3a (360_Cl, and 240_Ch). Figure 1 shows the HCV RNA levels and HCV-specific IgG antibody titers (estimated as the optical density to cut-off ratio in the enzyme immunoassay) over the course of infection. 10.1371/journal.ppat.1002243.g001 Figure 1 RNA level, Shannon entropy, and antibody titers over time for the four early infection subjects. Panels A and B show two subjects who developed chronic infection (240_Ch, 23_Ch) followed from pre-seroconversion timepoints. Panels C and D show the infection dynamics for two subjects who ultimately cleared the infection (360_Cl, 686_Cl). Red dots represent viremic time points analyzed via next generation sequencing (NGS). The solid line represents the RNA level. The dashed and dotted lines represent the interpolation of the Shannon entropy calculated across the genome at each time point using NGS data. Entropy was calculated using all mutated sites (dotted line), or with only non-synonymous sites (dashed line). The shaded area represents semi-quantitative estimates of the anti-HCV antibody titer (OD: cut-off). Note the varied ranges in the x- and y-axes. 10.1371/journal.ppat.1002243.t001 Table 1 Subject characteristics, number of founder viruses, and estimates of the time since the most recent common ancestor (tMRCA). Subject ID Sex Disease outcome GT Age First sampling timepoint (DPIa) Observed duration of viraemia (DPI) Timepoints analysed via NGSb Length of genome investigated Maximum entropy (full genomec) Number of founder viruses Poisson estimated tMRCA (95% CI) d 686_Cl F Cleared 1a 25 33 117 4 9172 0.010414 1 35 (26,47) 23_Ch M Chronic 1a 25 36 304 6 9138 0.002323 2 44 (24, 64) 360_Cl M Cleared 3a 29 30 223 2 5992 0.003190 1 34 (20, 57) 240_Ch M Chronic 3a 24 44 477 4 9226 0.001546 1 34 (19,48) a Estimated days post-infection (DPI). b Next Generation Sequencing. c Calculated across the genome for non-synonymous substitutions. d Mean value from the estimates of the time since most recent common ancestor (tMRCA) estimated from reconstructed viral variants across the genome in windows of 400 nt calculated according to a Poisson model of viral mutation (45). Confidence intervals are the 5th and 95th percentiles of the total range of confidence intervals estimated in each of the windows analyzed across the genome. For the purposes of this analysis, the time course of the primary infection was divided into three phases: i) transmission; ii) an acute phase (designated as 100 DPI) leading to clearance or chronicity at six months post-infection). Evolutionary dynamics of viral diversity over the course of the infection NGS, standard cloning and bulk (consensus) sequencing were performed on viremic samples collected longitudinally from each of the subjects (Figure 1). On average, 81,900 reads with average read lengths of 358 bp were generated per subject per timepoint, giving an average coverage depth of 3,093 (Table S2). Single nucleotide polymorphisms (SNP) analysis revealed between 160 and 460 nucleotide substitutions per timepoint. Over the course of the infection more than half of the substitution events occurred at frequencies of 99%). Colors represent the time course post-infection (see legend). Shannon entropy (SE) was calculated from the frequencies of nucleotide substitutions measured across the genome. During the acute phase, the SE measures indicated a general decline in viral diversity. This decline did not follow the HCV RNA kinetics. Similar patterns were observed for SE calculated with all substitutions, or with only non-synonymous changes. In two subjects (686_Cl, 23_Ch) at least two pre-seroconversion timepoints were sampled within a window of eight weeks, allowing an accurate assessment of very early HCV evolution, and revealing a peak in viral diversity in temporal proximity to the viremic peak. In the two subjects who ultimately became chronically infected (23_Ch and 240_Ch), SE increased in the pre-chronic phase, which also featured a continued increase in HCV RNA level (Figure 1). The observed patterns of change in SE measured across the full genome were not consistently observed in the 10 protein-encoding regions, clearly indicating a non-uniform evolution of HCV across the genome and between subjects. Transmission - First bottleneck To test the hypothesis that a single virus establishes infection upon transmission, a statistical model, PoissonFitter [44] was used to examine whether the viral population had a star-like phylogeny with a Poisson distribution [22]. For this analysis reconstructed viral variants were obtained from NGS reads using a Bayesian statistical tool, ShoRAH [34], [36]. Firstly, viral haplotypes were reconstructed from NGS data at the first viremic timepoint segregated into 400 nt windows across the genome (see Table S3, including the number of variants reconstructed per window). Secondly, reconstructed haplotypes of the E1/HVR1 (871 nt) and partial E2 (932 nt) regions were utilized, along with the corresponding E1/E2 clonal sequences (Table S3). In the second dataset, only variants with a frequency greater than 2.5% were included in this analysis. This cut-off was derived through a validation analysis for the method of reconstruction of viral haplotypes (of length > 400 nt) using a mixed sample of four plasmid E1/E2 clones derived from one subject (see Supporting Text S1 for details). By definition, the sequence of the founder was identified as being: i) the most prevalent variant; and ii) identical to the consensus sequence of the viral population at the first time point [15], [22]. Both the PoissonFitter tests and phylogenetic analyses indicated that the HCV infection in three subjects (686_Cl, 360_Cl, and 240_Ch) was successfully established by a single variant (Figures 3, S1, S2 and Table S3). For these subjects, the highlighter plots show the random distribution of SNPs across the sequences, again consistent with a star-like distribution of variants arising from a single founder. By contrast, Poisson analysis of the remaining subject (23_Ch) indicated that more than one virus established the infection. Phylogenetic analyses of E1, E2 and NS3 regions of subject 23_Ch indicated that at least two viruses had established the infection, designated 23AF and 23BF (Figures 3 and S1, Table S3). These founders generated two major clusters of variants, which were named after the founder variants. These two founder variants were identified with comparable prevalence in both E1/HVR1 (Figure S1) and E2 regions (Figures 3). The average genetic difference between these two dominant variants over E1/E2 regions was 1.3%. In support of these results, the founder variants were also identified via clonal sequencing in both E1/HVR1 (Figure S1) and E2 regions (Figure 3). For 23_Ch, the calculation of the exact number of founder viruses was confounded by the presence of several variants that appeared to have arisen from recombination events between two viruses from each cluster (see for example, variant C1_14b in Figure 3). 10.1371/journal.ppat.1002243.g003 Figure 3 Founder virus analysis based on partial E2 region of the viral genome. Panels A and B show the analyses for two subjects who developed chronic infection (240_Ch, 23_Ch) followed from pre-seroconversion timepoints. Panels C and D show the analysis for two subjects who cleared the infection (360_Cl, 686_Cl). Phylogenetic reconstructions and highlighter plots are shown, illustrating the genetic relatedness between HCV variant sequences. Names of each sequence are labeled with a letter (H for haplotype, and C for clone), with the first number representing the sampling timepoint and with the second number representing either the prevalence of the haplotype or the clone number. The phylogenetic trees of subjects 686_Cl, 360_Cl and 240_Ch (panels A, C, D) are consistent with an infection arising from a single founder. The fit with a Poisson model is also consistent with a single founder (p-value > 0.1, see text). As shown by the highlighter plots, founder viruses are identified as the consensus sequence and coincided with the most prevalent variant reconstructed from NGS data, (e.g. for subject 686_Cl H1_0.60 is identical to the consensus sequence and to clone C_12b). The highlighter plots also show the random distribution of mutated sites with respect to the founder sequence (master), which is consistent with a star-like phylogeny. The phylogenetic analysis in 23_Ch (panel B) is consistent with an infection originated from two founder viruses (indicated with an asterisk in the highlighter plot) giving rise to two major clusters, 23A and 23B. This is consistent with the rejection of the Poisson model (p-value = 0). Phylogenetic trees were obtained using PhyML, with Maximum Likelihood methods using a GTR model of substitution as suggested by model testing. A total of 2-5 variants (including the founders) were detected with a frequency above 2.5% in 240_Ch, 686_Cl, and 360_Cl at the first viremic timepoint. Subject 23_Ch presented the most diverse repertoire, with nine variants present in the E1/E2 region. These results, in combination with the fact that >50% of substitutions were present at a frequency 100 DPI) featured selective sweeps, with the new variants arising from the preceding viral population in the acute phase, showing remarkably few substitutions reaching fixation in the new population. In subject 240_Ch, following the decline in HCV RNA level at 140 DPI, a new cluster of variants, termed here 240AC, figure 4) replaced the 240AF cluster via a selective sweep brought about by the second bottleneck event. Similarly, in 23_Ch, the recrudescence of viremia in the pre-chronic phase was associated with the emergence of a new cluster, 23AC, which evolved only from cluster 23AF via a selective sweep (Figure 4b). In both subjects a single variant dominated the pre-chronic phase population, accounting for 74% of the total HCV RNA level in subject 23_Ch, and approximately 50% in 240_Ch. Further evidence in support of the establishment of a new viral population in the pre-chronic phase was provided by the increase in effective population size following the second bottleneck, which was particularly evident for 240_Ch (Figure 6). Fixation (>99%) was observed in a limited number of sites in association with the selective sweep (Figure 2). Fixation occurred at two sites in 23_Ch by 136 DPI (E2, 542; NS3, 1498) and at 12 sites in 240_Ch by 159 DPI (E1, 372; E2: 443, 483 and 543; NS2, 750; NS3: 1509, 1606; NS4b: 1768; NS5a: 1985, 1999; NS5B: 2497 and 2973). Further evolution was observed with additional fixation sites occurring by 304 DPI in NS2 (856) and NS5B (2629), indicating the presence of a further selective sweep event. In subject 240_Ch a further selective sweep event occurred by 477 DPI, with fixations observed in E2 (400, 405, 408, 583), NS3 (1641), and NS5a (2361). Evidence of immune driven selection In order to determine whether non-synonymous substitutions were likely to be driven by cellular and humoral immune responses, we identified the location of predicted, as well as experimentally-verified, epitopes across the full genome for each subject (Table S4). Analysis of HLA class I restricted cytotoxic T cell (CTL) epitopes revealed a total of 62 amino acid substitutions (29 in 23_Ch, 12 in 240_Ch, 8 in 686_Cl, and 13 in 360_Cl) which were found to be located in HLA-restricted epitopes (Table S4, Figure 5). Twenty-four of these substitutions resulted in mutant epitopes with a lowered predicted binding score, suggesting that these substitutions potentially conferred a CTL immune escape phenotype (Table S5). In 23_Ch, all four sites that reached fixation in the pre-chronic phase lay within HLA-restricted epitopes. In 240_Ch, of the 18 sites that reached fixation in the pre-chronic phase, nine lay within CTL epitopes. Figure 5 shows the location of the E2 mutations within putative CTL and B cell epitopes in the different variants in subjects 23_Ch and 240_Ch, as well as mutations previously shown to be associated with viral fitness costs. All the identified epitopes within this region carried at least one amino acid change. Two of these mutations (G483D for 240_Ch and T542I for 23_Ch) generated CTL epitopes with reduced binding affinity, thus indicating potential immune escape (Table S4). Interestingly, both subjects also showed a substitution at position 443 (Y443I for 23_Ch and Y443H in 240_Ch) known to be within a B cell epitope [47]. In 240_Ch a mutation at 543 was observed (T543A). Mutations at this site have been reported to alter the yield of precipitated E1/E2 proteins [48]. Comparison between SE (based on all substitutions) at sites of predicted CTL epitopes, with sites at which no epitopes were predicted, revealed significant differences for each subject (Mann Whitney; all p 2.5% were considered for further analysis (see Supporting Text S1). The probabilistic nature of the clustering algorithm allows for estimation of the reliability of predictions. In the Bayesian scenario implemented in ShoRAH, the fraction of iterations in which a haplotype is reported estimates the posterior probability of the existence of that haplotype. This posterior probability provides a confidence level for the haplotype. Only haplotypes with confidence values (posterior probabilities) greater than 0.9 are reported. Scripts and parameter values used for the analyses are available upon request. Parameters of ShoRAH have been chosen through a detailed preliminary analysis on the sensitivity of these parameters in collaboration with Dr Zagordi (personal communication). Founder virus analysis PoissonFitter was used to test the hypothesis that a single virus establishes infection [71]. PoissonFitter performs two tests: one test is based on the fit of the Poisson model to the frequency distribution of the Hamming distance observed in each sample; the other is a topological test to verify that observed frequencies are distributed according to a star-like phylogeny (for this test, no formal statistic is available and consequently no p-value is obtained). In this model the main assumption is that a single founder virus evolves under neutral evolution, generating a star-like phylogeny, with a distribution of mutations conforming to a Poisson distribution [15], [45]. This means that early selective pressures compromise the statistical analysis. For this reason, only HCV sequences obtained from the first available viremic time point for each subject were included in the founder virus analysis. PoissonFitter test was performed on reconstructed viral variants in E1/E2 region of the genome, as well as on reconstructed variants obtained from 400 nt windows sliding across the full genome. A type 1 error threshold α = 0.01 was assumed to conclude whether the Poisson test was rejected. For the initial analysis, in each window a conclusion was made whether one, or more than one, founder virus explained the observations. In cases where the two tests described above were discordant (i.e. a significant fit with Poisson test of a single founder virus, but no star-like phylogeny – which occurred in seven analyses; or the Poisson test suggested more than one founder while the topological test conformed to a star-like phylogeny – which occurred in eight analyses) then early stochastic events were assumed to have occurred, which limited the validity of the tests. These early events could result from an early selective pressure, or strong fitness advantage associated with early mutations. In this analysis, the mutation rate for the HCV genome was set to 1.2×10−4, as recently estimated [72]. We assumed the same value was constant across the HCV genome. Phylogenetic analysis Sequences from standard cloning and reconstructed haplotypes were visualized and curated with MEGA [73] and R packages. Phylogenetic and evolutionary analyses were performed with PhyML [74]. Trees were constructed from sequences using a GTR substitution model with gamma invariant sites, as suggested by analyses using ModelTEST [75]. Trees were visualized with FigTree. Estimates of tMRCA were performed with PoissonFitter and with BEAST [76]. Demographic reconstruction of the within-subject viral population was performed with skyline plots implemented in BEAST [47], [63]. The analysis was performed using both a strict clock model and a relaxed molecular clock model. The length of the MCMC chain was chosen so that the effective sample size (ESS) for each parameter was > 100. Bayesian factors were then used to decide the best model. For each subject, the null hypothesis of a strict clock was not rejected according to the Bayes Factor calculated from the posterior distributions obtained from each model (BF 0.1, see text). As shown by the highlighter plots, founder viruses are identified as the consensus sequence and coincided with the most prevalent variant reconstructed from NGS data, (e.g. for subject 686_Cl H1_0.70 is identical to the consensus sequence and to six clone sequences, see Highlighter plot). The highlighter plots also show the random distribution of mutated sites with respect to the founder sequence (master), which is consistent with a star-like phylogeny. The phylogenetic analysis in 23_Ch (panel B) is consistent with an infection originated from two founder viruses (indicated with an asterisk in the highlighter plot) giving rise to two major clusters, 23AF and 23BF. This is consistent with the rejection of the Poisson model (p-value = 0). Phylogenetic trees were obtained using PhyML, with Maximum Likelihood methods using a GTR model of substitution as suggested by model testing. (TIF) Click here for additional data file. Figure S2 Fit of the distribution of Hamming distances with a Poisson model performed on the Envelope region, specifically on partial E2 (932 nt). Each panel shows a subject; the red line is the fit of the distribution of Hamming distances (histogram) calculated from sequences derived from the reconstruction of NGS data and from standard cloning. The poor fit outcome for subject 23_Ch indicated the presence of more than one founder virus. See also Table S3. (TIF) Click here for additional data file. Figure S3 Evolutionary dynamics of HCV variants over the E1/HVR1 region of the genome. Sequence analyses of the two subjects who developed chronic infection, 240_Ch (A), and 23_Ch (B) revealed the presence of selective sweeps. These sweeps led to the emergence of new variants that replaced the founder viruses (identified with an asterisk). Phylogenetic trees (left panels) display nucleotide sequences of reconstructed haplotypes derived from NGS data and clonal sequences. Names of each sequence are labeled with a letter (H for haplotype and C for clone), with the first number representing the sampling timepoint and with the second number representing either the prevalence of the haplotype or the clone number. Colors are also used to portray the sampling timepoint (see legend). Infection dynamics for subject 240_Ch are consistent with a single founder, identified with the most prevalent strain of cluster 240A (H2_0.65, indicated with an asterisk), and with clones C2_6, and C2_17, and with the consensus of the sequences from time-point 1. The pre-chronic phase (corresponding with the color-coded time ranges 5 and 6) of infection shows the emergence and dominance of a new subgroup of viruses. 23_Ch has two founder viruses that successfully initiated the infection (H1_0.25 and H1_0.27, indicated with an asterisk). A new cluster emerged from the founder cluster and became dominant in the pre-chronic phase. Trees are calculated using Maximum Likelihood method (implemented in PhyML). (TIF) Click here for additional data file. Figure S4 Evolutionary dynamics of HCV variants over the partial NS3 region of the genome. Sequence analyses (over 400 nt window) of the two subjects who developed chronic infection, 240_Ch (panel A), and 23_Ch (panel B) revealed the presence of selective sweeps. These sweeps led to the emergence of new variants that replaced the founder virus(es). Phylogenetic trees display nucleotide sequences of reconstructed haplotypes derived from NGS data and clonal sequences. Names of each sequence are labeled with a letter (H for haplotype), with the first number representing the sampling timepoint and with the second number representing the prevalence of the haplotype. Colors are also used to portray the sampling timepoint (see legend). (TIF) Click here for additional data file. Figure S5 Demographic reconstruction of the viral populations. Demographic reconstruction from E1-HVR1 and E2 sequences for subjects who cleared the infection (686_Cl and 360_Cl). In both subjects and in both genomic regions the estimated effective population size (Nτ, the product of the effective population size and generation length in days) is higher than in chronic subjects with peak around 104. Subject 686_Cl shows a moderate increase in viral diversity over time, corroborating resuts from Shannon entropy and Poisson-Fitter tests. (TIF) Click here for additional data file. Table S1 Demographic and laboratory characteristics of the subjects. (DOC) Click here for additional data file. Table S2 Summary of the next generation sequencing (NGS) data obtained for each subject. (DOC) Click here for additional data file. Table S3 PoissonFitter test results for the single founder virus analysis for each subject. (DOC) Click here for additional data file. Table S4 Phenotypic analyses of the non-synonymous substitutions over the course of the infection. (XLS) Click here for additional data file. Table S5 Analysis of amino acid substitutions within predicted HLA restricted cytotoxic T cell epitopes. (DOC) Click here for additional data file. Table S6 Primers used to amplify viral sequences. (DOC) Click here for additional data file. Text S1 Haplotype reconstruction and method validation. Detailed information on the haplotype reconstruction from NGS data. This methodology was validated with a 454 FLX run performed in duplicate on four unique plasmids with known sites of variation. (DOC) Click here for additional data file.
