Antigenic drift and the generation of viral quasispecies

Some RNA viruses form a quasispecies--a set of related viral variants that coexist in field populations and even within single infected individuals (reviewed in [1,2,3,4,5]). The emergence of immunologically distinct members of a viral quasispecies through mutation and subsequent immune selection is called "antigenic drift." Antigenic drift is thought to be important in human immunodeficiency virus (HIV) infection and the continuing seasonal influenza epidemics because immunity generated against one viral quasispecies member selects for escape variants. Attributed in part to antigenic drift are the moderately high failure rate and the short-lived efficacy of influenza vaccines [6], the failure of synthetic foot-and-mouth disease virus vaccines [7], and the inability of recombinant HIV vaccines to provide complete protection against field strains of the virus [8].

The hemagglutinin (HA) envelope surface glycoprotein--the major neutralizing determinant of influenza A--is a classic example of an antigenically drifting protein [9]. Walter Gerhard and colleagues demonstrated that the immune pressure exerted by monoclonal antibodies (Abs) selects for HA escape mutants in model systems [10,11]. Later, Dimmock and colleagues showed that polyclonal anti-sera also select for escape mutants [12,13]. Similarly, much of the observed variability of glycoprotein 120 (gp120), the principal surface antigen of HIV, is thought to reflect antigenic drift [14,15,16,17]. The correlation of intra-patient viral diversity with immune response strength has been cited as evidence that the immune response is a selective factor in HIV antigenic drift [18,19,20,21,22].

Phylogenetic analyses describe divergence within a viral population, and these methods have been used to infer the selective advantages of viral variation [18,19,20,21,22,23,24,25,26,27,28,29,30]. A more direct indication of the selective advantage gained through variation is an observed overabundance of replacement mutations relative to silent mutations in viral proteins [31]. Such analyses of gp120 and its V regions indicate that replacement mutations are generally over-represented in this protein and thus appear to confer selective advantage to HIV-1 [22,24,25,26,27,28,29,30,32,33,34,35]. In more detailed analyses, several groups tested individual codons for replacement mutations that are, as an aggregate, overabundant [23,36,37,38]. However, none of these methods determine which replacement mutations are actually positively selected. Also, when replacement mutations of varying fitness are lumped together, positively selected mutations may remain undetected among negatively selected mutations.

To overcome these limitations, we have developed a "selection mapping" algorithm. The cornerstone of selection mapping is the testing of each observed replacement mutation at each codon to identify those particular replacement mutations that are overabundant relative to silent mutations at that codon. Such replacement mutations are determined to be positively selected. Negatively selected variants are recognized as "noise" and are thereafter ignored. Here, we use the selection mapping method to identify the positively selected variants of influenza A HA (H3 serotype), HIV-1 reverse transcriptase (RT), and HIV-1 gp120.

Selection map of influenza A H3 hemagglutinin

QUASI identifies 25 HA codons where one or more replacement mutations are positively selected in the influenza A H3 virus (Fig. 1). [From our neutral drift testing of QUASI, we expect a maximum of about 2-3 false positives (see Materials and Methods)]. The distribution of these positively selected codons is of particular interest. Without exception, the codons where variants are positively selected are on the HA surface (Fig. 2, left). The parsimonious explanation of this result is that HA variants are primarily selected for escape from B-cell immunity. If T-cell immunity is instead the primary selective force affecting HA variation, then either all T-cell immunity escape variants are coincidentally solvent-exposed, or HA T-cell epitopes are determined by Ab protection [39]. These 25 codons include 13 outside those identified as positively selected by Walter Fitch and colleagues [40]. We attribute our new findings primarily to sites where many variants are negatively selected but where at least one variant is positively selected. The Fitch group also identifies 6 additional codons as positively selected. We believe that these are false positives caused by the Fitch group's assumption that HA is, on average, neutrally drifting; these six codons may be less negatively selected than average, but they nevertheless appear to be negatively selected. Additionally, our data, while largely the same as those analyzed by the Fitch group, do have some differences.

Wiley et al. proposed four antigenic sites where field and laboratory mutations could be grouped on the HA surface [41]. These putative antigenic sites are indicated in Figure 1 and the right side of Figure 2. Positively selected variants are correlated (Fisher's exact test) with antigenic site A (p = 5.27 × 10^-3), antigenic site B (p = 1.12 × 10^-3), and antigenic site C (p = 1.0 × 10^-5). In contrast, antigenic site D is not particularly correlated with positively selected variants. We believe the lack of positively selected variants spanning antigenic site D may explain a decades-old puzzle. In the fully assembled HA protein, site D is buried in the trimer interface and therefore is not generally accessible to Ab [41,42]. At the time that the antigenic sites were proposed, residues in and around the trimer interface were variable and found on the monomer surface, so they were grouped and labeled as antigenic site D even though it was unclear how site D was recognized by Ab [41]. In fact, the only two positively selected variants QUASI identifies in site D are solvent-exposed at the extreme edge of the HA trimer (Fig. 2). Based on QUASI's selection map of HA, we now conclude that those mutations found in the trimer-buried portion of HA do not confer significant advantage on influenza A. That is, while the trimer-buried portion of antigenic site D includes the sound and fury of variability, it signifies nothing.

While QUASI finds that antigenic sites A-C are significantly associated with positively selected variation, this association may be a simple consequence of the sites' surface exposure. As can be seen in the left-hand side of Figure 2, positively selected variants are scattered across the entire exposed surface of HA. There are positively selected variants outside the antigenic sites, and there are subregions of the antigenic sites where variation is negatively selected. Thus, it may be more appropriate to view antigenic sites as more positional than functional. Indeed, more recent demarcations of antigenic sites enlarge antigenic sites A-D and add an additional antigenic site, site E, and these sites are now considered primarily positional rather than functional [6,43].

Selection map of HIV-1 gp120

At 123 codons, QUASI's selection map of gp120 indicates that one or more non-consensus amino acids are positively selected in HIV-1 (Fig. 3). Most positively selected variants appear to be on the gp120 surface (not shown), but in contrast to HA, gp120 includes several monomer-buried positively selected variants (at sites I225, V270, N295, H333, I345, T387, I424, and L453). Additionally, some of the positively selected variants that appear to be solvent-exposed may normally be buried (some loops are absent from the core protein crystal structure [44]). The burial of positively selected variants in the gp120 monomer confirms that gp120 quasispeciation is not selected solely for escape from B-cell immunity.

Two competition-group epitopes have been identified for broadly neutralizing anti-gp120 Abs: the CD4-binding site (CD4BS) and the CD4-induced (CD4i) epitopes (references in [45]). Each epitope includes only a single non-consensus positively selected variant (Fig. 3). Thus, broadly neutralizing Abs appear to be those that engage few protein positions where variation is positively selected. Based on this observation, we propose that the neutralizing spectra of Abs may be predicted if the epitopes are known. For example, we would predict that the anti-CD4BS Abs will have particular trouble recognizing gp120 molecules carrying the positively selected variant of the CD4BS epitope (D→N at codon 474). This prediction is fulfilled in that the 15e anti-CD4BS monoclonal Ab fails to react to gp120 from HIV strain RF, a strain that carries this D→N mutation [46]. We also predict that gp120 molecules carrying the positively selected I→F mutation at codon 423 will be poorly recognized by the anti-CD4i Abs.

It is worth commenting on the GPGRAF motif (gp120 residues 312-317) that is sometimes (though increasingly rarely) referred to as "highly conserved." Because QUASI identifies non-consensus variants at two codons of this motif as positively selected (Fig. 3), it may be inappropriate to refer to the GPGRAF motif as "highly conserved." Instead, five positively selected variants appear to exist at this region in addition to GPGRAF: GPGKAF, GPGRTF, GPGKTF, GPGRVF, and GPGKVF.

Kwong et al. roughly divide the gp120 three-dimensional structure into outer (β9-β19 and β22-β24) and inner (N-α1, β4-β8, and α5-C) domains joined by a bridging sheet (β2, β3, β20, and β21) [44]. As indicated by Kwong et al., all three domains include variable regions; as QUASI shows, diversity in all three regions is positively selected (Fig. 3). The selective advantage we find rendered by diversity in some of these regions (e.g., the "silent face") has been attributed to neutral drift [44,45], thus QUASI's results run counter to previous interpretations. That is, QUASI finds that diversity in regions proposed to be inaccessible to the immune system nevertheless confers selective advantage on HIV. QUASI's findings may be consistent with the existence of gaps in the carbohydrate groups thought to mask the silent face from gp120 from immune surveillance. Interestingly, in carbohydrate-building models of gp120, these gaps correspond to codons where we find variation is positively selected (P. Kwong, pers. commun.).

A "non-neutralizing" face has been identified where binding Abs generally do not neutralize HIV when gp120 is oligomerized [45,47]; these data were interpreted as indicating that the non-neutralizing face is occluded in the trimer and that binding Abs are raised against shed gp120 monomers. However, QUASI finds numerous positively selected variants on the non-neutralizing face of the inner domain. Assuming the "non-neutralizing" appellation is appropriate, how do mutations on this face provide selective advantage to the virus? The obvious answer is that mutations provide escape not from direct B-cell immunity but from other levels of immunity, such as T-cell immunity [including major histocompatibility complex (MHC) presentation] or indirect Ab immunity via Ab-dependent cellular cytotoxicity (where gp120 molecules found on infected cell surfaces are monomers).

To determine if HIV-1 viral sequences retain evidence that T-cell immunity is a significant selective force affecting HIV quasispeciation, we used QUASI to generate a selection map of HIV-1 RT (Fig. 4). Because RT is not a surface-expressed protein, it is not plausible that the positively selected variants of RT have been selected by direct B-cell immunity. A priori, RT quasispeciation could have been the result of neutral drift, but because QUASI finds that replacement mutations confer selective advantage on the virus, we reject the neutral drift hypothesis at the 22 RT codons. If T-cell immunity (including MHC presentation) is a selective pressure shaping RT quasispeciation, positively selected variants should be associated with T-cell epitopes. When known T-cell epitopes [48] are plotted on the RT selection map, the positively selected variants are found to localize significantly (Fisher's exact test) both with helper T-cell epitopes (p = 3.27 × 10^-2) and CTL epitopes (p = 6.58 × 10^-3). We conclude that T-cell immunity is a significant selection pressure shaping the quasispeciation of RT and presumably is a significant factor in the quasispeciation of other HIV proteins.

Thus, because positively selected gp120 variants found throughout gp120 may be selected by T-cell immunity, QUASI's finding that the non-neutralizing face includes positively selected variants is not at odds with models where the non-neutralizing face forms the gp120 trimer interface. Nor are QUASI's results incompatible with the silent face being silent to B-cell immunity. QUASI's finding that positively selected variants may be buried in the gp120 monomer is consistent with escape from T-cell immunity.

In addition to the selection pressure exerted by T-cell immunity, 3'-azido-3'-deoxythymidine (AZT) may also have provided selection pressure for RT quasispeciation in the sequences selection mapped by QUASI [59,51]. Indeed, QUASI identifies six of the eight mutations known to be associated with AZT resistance [52] as positively selected (N67, R70, W210, Y215, and F215) or possibly positively selected (L41) to HIV-1 (Fig. 4). The exceptions, two mutations of codon 219, are informative. Whereas mutations at other codons are necessary for high resistance to AZT, mutations at codon 219 are not, and codon 219 mutations arise late in infection after earlier mutations have already rendered RT resistant to AZT [53]. We conclude that the additional AZT resistance conferred by codon 219 mutations did not provide significant selective advantage to the profiled HIV viruses, possibly because HIV had already acquired the maximum effective AZT resistance selectable, in vivo, when mutations arose at this codon. Alternatively, the lysine at codon 219 may be important for proper in vivo RT function such that the advantage conferred by increased AZT resistance does not adequately compensate for impaired RT function. The RT sequences we analyze were taken from patients who either had no anti-RT treatment or were treated mainly with AZT (though some patients who were treated with AZT were also treated with 2',3'-dideoxyinosine) [49,50,51]. Therefore, we would predict that mutations associated with resistance to other anti-RT drugs should not be positively selected in the sequences QUASI analyzed. As predicted, QUASI identifies none of the 50 RT mutation associated with resistance to other drugs as positively selected (compare Figure 4 to the Los Alamos database [52]).

QUASI--the selection mapping algorithm

An executable version of the QUASI software is attached as an additional file (see additional file 1). Also attached are a users' manual (user.txt - see additional file 2) and a FASTA to QUASI file converter PERL script (F2Q.pl - additional file 3). Current versions of QUASI are available from the authors or may be accessed at the Los Alamos Influenza Sequence Database (http://www.flu.lanl.gov/).

For a set of viral nucleotide sequences, we determine the variants that confer selective advantage by measuring the empirical replacement to silent mutation ratio (R:S) of each possible amino acid replacement and then comparing this observed ratio to that which would be expected if mutation were unselected. An R:S that is found to be higher than expected indicates that the replacement mutation tested is positively selected, while a lower-than-expected observed R:S indicates that the tested replacement mutation is negatively selected.

Testing for an overabundance of replacements across a protein as a whole is a reasonable approach when only a few nucleotide sequences are available, but because a large number of mutated viral sequences are currently available, such aggregation is unnecessarily crude. Better are approaches that test for an overabundance of replacement mutations at individual codons [23,36,37,38]. However, these methods lump together replacement mutations and thus allow negatively selected mutations to conceal positively selected mutations and vice versa (e.g., replacement mutations at a codon may be negatively selected as a group despite the fact that one or more particular replacement mutations are positively selected).

To overcome these limitations, the QUASI algorithm does not test the overall R:S of the entire protein as an aggregate, nor does QUASI test the R:S of a codon to all its replacement mutations taken as a whole. Rather, QUASI tests the R:S of each particular replacement mutation at each codon. That is, QUASI measures the R:S of the mutations from a consensus codon towards each individual replacement amino. For example, if the consensus codon at a protein position were ttt (Phe), QUASI would test the R:Ss of all point mutations from ttt. One of these mutations is ttt→tat (Tyr). QUASI calculates the expected R:S for ttt→tat under the null hypothesis of neutral drift. The expected S is one because only one mutation (ttc) is silent, and the expected R is also one because only one point mutation of ttt (tat) codes for Tyr, so the expected R:S is one, in this case. If QUASI rejects the (Jukes-Cantor) neutral drift null hypothesis because the observed R:S is significantly higher than one, then QUASI classifies this replacement mutation (Tyr) as positively selected. Conversely, if QUASI rejects the null hypothesis because the observed R:S is significantly lower than one, then QUASI determines that this replacement mutation is negatively selected. QUASI performs this procedure for all replacement point mutations [e.g., in the example case, Tyr (tat), Ile (att), Leu (tta, ttg, and ctt), Val (gtt), Ser (tct), and Cys (tgt)].

In this paper, selection mapping is carried out independent of the underlying phylogeny. QUASI uses R:S to reject the null hypothesis that the mutational space surrounding the consensus codon is distributed randomly among all nine possible R or S point mutations (except stop codons, which are considered to be disallowed). This allows R:S calculations to be applied to viral sequences whose ancestral sequence is unclear or unknown. This is a both an advantage and a disadvantage over analyses that rely on phylogeny. Phylogeny is difficult to determine accurately and uniquely, and relying on phylogeny ignores the persistence of positively selected replacement mutations (the major effect of selection). On a practical level, using phylogeny to reconstruct viruses' mutational histories and then using intuited mutations leaves one with insufficient data to determine positively-selected codons [36] unless, as some have done, one assumes observed drift is neutral and then tests for codons where selection is more positive than average [23,38]. The significance problem can be compounded when one is looking for independent occurrences of particular mutations; often, there simply has not been enough sequence evolution in HIV or influenza to map positively selected variants if the retention of positively selected mutations is ignored. The drawback of ignoring phylogeny is a potentially high false positive rate (see below).

Empirical R:S is compared to neutral R:S by means of a two-sided test of the binomial distribution. For each codon, we test the null hypothesis that all nine point mutants are equally probable. The quotient p = R/(R+S) is the probability of a replacement mutation at this codon if each nucleotide is equally mutable and each of the three mutational targets at that codon are equally likely. The numerator, R, is the number of point mutations that lead from the consensus codon to the target amino acid. The chance of observing r replacement mutations is given by the binomial distribution, , where n is the number of codons providing data for this position. To form a two-sided test, we sum all terms b(kn,,p) such that b(kn,,p) is not greater than b(rn, p,), wherek is in the set (0,...,n) and r is the number of observed replacement mutations. In other words, we sum the chances of all events that are no more likely than that of the observation. If this sum, α, is small (e.g., not greater than 0.05), we reject the null hypothesis at the α level of significance.

Working example

The consensus codon is given as ttt (Phe; Table 1, column 1). Each observed mutation is also given (Table 1, columns 1-3)

Because we know the frequency of silent mutations (given as 10 in this example; Table 1, column 3), we also know the expected R:S for each replacement mutation (Table 1, column 4). That is, if selection is neutral for any particular replacement mutation, we can calculate the incidence of each replacement mutation we expect to observe (by looking at a table of the genetic code).

Using the given frequency of each mutation observed, we also know what the observed R:S is in each case (Table 1, column 5).

Now we use a two-tailed test of the binomial distribution to determine if each observed R:S is significantly different from the corresponding expected R:S (Table 1, column 6). In some cases, the differences are significant, in which case the appropriate replacement mutation is assigned positive or negative selection (positive if the observed R:S is larger than expected and negative if the observed R:S is lower than expected). Otherwise, the selection is assigned to be neutral drift.

The QUASI algorithm thus indicates at this exemplary codon that both tyrosine and serine are positively selected; isoleucine, leucine, and cystine are negatively selected; and the selective advantage or disadvantage of valine is indistinguishable from neutral drift. Any other ttt codon will have its own selective pressures assessed in a similar but independent testing procedure.

Minimum sequences

For each possible replacement mutation, a minimum number of mutations will need to be observed for a selective event to be detected. This minimum number differs depending on the consensus codon and the level of significance. We have calculated the minimum number of mutations needed to achieve the 5% significance level. At the lower bound of this range, only 2 replacement mutations will be required to detect positive selection for any replacement mutation from cta or ctg. Any replacement mutation from cta or ctg has an expected R:S of 1:4, and thus an observed R:S of 2:0 will be sufficient to reject neutral drift in favor of positive selection. Conversely, detecting negative selection at a cta or ctg codon is difficult. At the upper bound of the range, a minimum of 17 observed mutations are required to detect negative selection at such a codon (0:17 is significantly lower than 1:4). For the most-typical codon, the expected R:S is 1:3. For these modal codons, at least 3 mutations must be observed to detect positive selection (i.e., if observed R:S = 3:0). At the same most-typical codon, 12 mutations are required, at a minimum, to detect negative selection (i.e., if observed R:S = 0:12). Because identification of positive selection is generally the goal, the QUASI algorithm appears to have a practical advantage over extant selection detectors, which require either many more mutations or a biased expectation of neutral drift in order to detect positive selection.

False-positive testing

False positives are likeliest when drift is completely neutral (e.g., as was found by Suzuki and Gojobori [36]). One may estimate the maximum frequency of false positives by testing simulated sequences generated under neutral drift parameters. We used the EVOLVER program of the PAML package [54] to generate sequences drifting under neutral Jukes-Cantor evolution. For each simulation, we generated 300 related sequences of length 999 and with average branch lengths varying in 0.1 length intervals from 0.1 to 1.0; each parameter set was used to generate 10 sets of 300 sequences. False positive percentages [Fig. 5; false positive percentage = false positives / (false positives + neutral drift variants)] were extremely low (∼ 2%) for relatively long branch lengths (0.1-1.0). Extremely high false positive rates (up to 70%) were found for extremely short branch lengths (maximum at 0.001). As accurate branch lengths are calculated using maximum likelihood phylogeny, these branch lengths may be used to find the appropriate false-positive percentage. For instance, if the branch length were 0.01 [as appears appropriate for HA (unpublished observation)], then the maximum false positive rate would be estimated at 32%. We then use the calculated neutral drift frequency (from QUASI) as an estimator for the maximum false positives. If (as we report) HA is found to have 5 neutral drift variants, then we would expect a maximum of 2.3 false-positives.

Selection mapping

Sequences

Nucleotide sequences were downloaded from GenBank at the NIH. Sequences were included only if they were isolated in the field and were not obviously pseudogenes (sequences containing premature stop codons were removed from consideration). We analyzed 310 sequences of human-infective influenza A (H3 serotype), 6,151 HIV-1 gp120 sequences, and 400 HIV-1 RT sequences. All sequences were pre-aligned with PILEUP [56] and/or DIALIGN2 [57] then hand-corrected.