Determination of BL.

A maximum likelihood (ML) tree was constructed previously for each transmission pair (18), with gapped regions (including sequences in the first, second, and fourth hypervariable domains [V1V2 and V4]) stripped out of the analysis. Modeltest was used to identify the optimal evolutionary model (60), and for each tree, the model selected included base frequencies, a general reversible model for base substitution, and differences in rate variation at different sites assigned according to a gamma distribution with four categories of rates and invariant sites. A further comparison was made with a REV model with ML assigned rate variation at different sites, but the additional parameters were not justified in a likelihood ratio test for these data, and thus the Modeltest model was used. The specific value for each of the model parameters was estimated for each of the eight trees separately. PAUP was then used to construct a likelihood tree for each transmission pair. BranchLength.pl (B. Korber [www.santafe.edu/~btk/science-paper/bette.html]; a web-based version of this code is now available at http://hcv.lanl.gov/content/hcv-db/BRANCHLEN/branchlength.html) was used to calculate the branch length (BL) to the patient's ancestral node for each donor and recipient sequence.

Correlation analysis.

A nonparametric Kendall tau rank correlation test was applied to previously described paired sequence and neutralization data (18) to determine whether the 50% inhibitory concentration (IC₅₀) of plasma was correlated with the following different parameters: BL to the patient's ancestral node on the ML tree, amino acid length of the V1-to-V4 region, and the number of N-Gly sites in V1-to-V4. P values of <0.05 were considered statistically interesting for hypothesis forming as a foundation for further experimental exploration. The correlation tests were performed with the donor plasma IC₅₀ and pooled plasma IC₅₀ against each sequence attribute.

Mutual information analysis.

Shannon entropy is a measure of uncertainty used in communication theory (64). It can be applied to other problems where character states are considered, such as sequence analysis (37). For example, for position i in a protein alignment, the Shannon entropy function H(i) is defined by the probability of different amino acids, s, in that position, as follows: H(i) = −S·P(s_i)·log P(s_i). This is a measure of the uncertainty in the column of data; if only one amino acid appears in a column, it has zero entropy. Maximal entropy occurs when the frequencies of all of the characters are equal. Here we considered j a character state assigned according to the level of neutralization sensitivity associated with a V1-to-V4 sequence (a, IC₅₀ of <50; b, IC₅₀ of 50 to 249; c, IC₅₀ of 250 to 999; d, IC₅₀ of 1,000 to 4,999; and e, IC₅₀ of >5,000). We asked whether j covaries with amino acid patterns found in any given position in an alignment. Given the joint probability distribution P(s_i,s_j), H(i,j) = −S·P(s_i,s_j)·log P(s_i,s_j) and the mutual information distribution M(i,j) = H(i) + H(j) − H(i,j) (37), we measured how much the data covaries, given that the mutual information is 0 when either i and j are completely independent or there is no variation. To test the significance of these results, 10,000 randomizations of neutralization phenotype assignments to viral sequences were performed, and the number of times the mutual information for the random data was greater than or equal to the original data was tallied and divided by 10,000 to estimate the P value. This statistical approach is based on the assumption that the data are drawn from a random pool; however, the sequence data have strong underlying phylogenetic correlations (2). To partially compensate for this, the randomization of neutralization assignments to sequences was done within donor-recipient pairs; however, the associations identified must therefore still be considered with this limitation in mind and thus viewed in a hypothesis-forming framework to guide experiments and for comparison with other results, such as the dN/dS ratios. Hypervariable loops were omitted from this analysis because we could not exclude the possibility that alignment artifacts could occur in the association patterns.

Codon-specific dN/dS ratios.

Previously described data regarding codon-specific dN/dS ratios were adapted (22). Nonsynonymous mutations become fixed with greater or lesser probability than synonymous mutations depending on whether a site is subject to positive or purifying selection, and thus the dN/dS ratio can be used to characterize specific codons and regions of proteins evolving under positive selection (77). In viruses, regions under strong positive selection occur where immune escape confers a selective advantage to new variants that arise within a host. Accordingly, codon-based models of molecular evolution can be used in conjunction with sequence data and ML parameter estimation to identify sites of potential immunological interest. The dN/dS ratio is indicative of the selection pressure at the protein level, as follows: a dN/dS ratio of <1 is indicative of purifying selection and amino acid conservation because of structural and functional constraints, and a dN/dS ratio of >1 is indicative of diversifying, positive selection where amino acid substitutions confer an advantage. In this study, we asked whether or not the precise positions under positive selective pressure may be lineage specific for HIV-1, comparing subtype B and subtype C, as representatives of both lineages are under consideration as vaccine candidates in trial populations where the HIV-1 C clade dominates. The program CODEML (S8 [http://abacus.gene.ucl.ac.uk/software/paml.html]) was used to determine the regional distributions of dN/dS ratios in the V3 loop and flanking regions for two sets of 25 HIV-1 sequences, one for subtype B and the second for subtype C. The subtype C sequences used did not include those from the Zambian transmission pairs analyzed for neutralization sensitivity here. For each subtype, a three-rate prior with flexible locations and weights was fitted to the data. Bayes's formula was then used to calculate the posterior distribution of the rate ratio at each codon. The posterior mean was taken as a site-specific estimate of the dN/dS ratio.

Three-dimensional mapping.

The core structure used for mapping corresponds to the X-ray structure of CD4-bound YU2 gp120 (43) (PDB code 1RZK). Variable loops V1V2 and V3 were modeled for clarity as described previously (5). Discussions on an unliganded gp120 core structure were based on the fully glycosylated simian immunodeficiency virus (SIV) structure (12) (PDB code 2BF1). Signature positions were mapped onto this structure based on the alignment of sequences. Modeling calculations were performed using a modified version of AMBER (56). Three-dimensional images were generated using VMD (29).

Molecular dynamics simulations.

Comparative all-atom molecular dynamics simulations were carried out on the modeled gp120 to capture the local differences in the V3-C3-V4 regions of gp120 between B and C subtypes. The following three sets of simulations were carried out under identical conditions: (i) gp120 from YU2 (subtype B), (ii) YU2 gp120 with the α2 helix consensus sequence from subtype C, and (iii) gp120 with the subtype C consensus sequence. In all cases, the gp120 protein was solvated in ~16,000 water molecules. The starting structure for subtype C gp120 was based on the YU2 model described above. The subtype C consensus was homology modeled from the subtype B structure, using MODELLER software (69). All molecular dynamics simulations were performed using the force field described by Cornell et al. (15) and the AMBER 6.0 suite of programs (9). A rectangular box of TIP3P water molecules (31) was added to solvate the complex, with the final system containing ~55,000 atoms within a box dimension of ~79 by 79 by 89 Å. Periodic boundary conditions were applied, and the particle mesh Ewald approach (16) was used to accurately treat electrostatic interactions. Both systems were minimized and equilibrated with the same protocol. The temperature was maintained at 300 K, using the Berendsen temperature algorithm (3). A 2-fs time step was used, and all bonds involving hydrogen atoms were constrained using the SHAKE algorithm (68). Simulations were carried out for approximately 4 ns, and the last 2 ns were considered for analysis. Additional details of the simulations and the results are provided in a separate publication (22a). The relative solvation between subtypes was obtained by calculating the number of water molecules in the first solvent shell surrounding a secondary structure (e.g., the α2 helix) or a specific residue. The contact distances between C-alpha carbons were used to identify interactions between regions or differences between the clades. A 10-Å cutoff between C-alpha carbons was used as a measure for contact.

Construction of chimeric Envs.

PCR amplification and cloning of the donor and recipient Envs from pairs 55 and 135 have been described previously (18). The env genes are in the cytomegalovirus-driven expression vector pCR3.1 (Invitrogen), which is used to generate viral pseudotypes. Each of the donor α2 helix domains shown in Fig. ​Fig.3A3A was inserted into the corresponding recipient Env by using a strategy similar to that described previously (66). A DNA duplex corresponding to each donor α2 helix was synthesized as a double-stranded oligonucleotide by IDT Technologies (Coralville, IA) and gel purified to ensure that each fragment was of the correct length. These fragments were each blunt end ligated to an ~8-kb PCR amplicon corresponding to the recipient Env, excluding the α2 helix and including the pCR3.1 vector sequences. The single-stranded DNA sequences corresponding to the 135 F and 55 M donor α2 helices (HXB2 nucleotides [nt] 7227 to 7279) were as follows: for pair 135, F42a (5′-AAAAAGCTATGGAATAACACTTTAGACGGGATAAAGGAAAAATTAGCAAGATAC-3′), F48a (5′-AAACAAAATTGGACTAAAACTTTAGAAGAGGTAAAGACAGAATTAAGAAAATAC-3′), and F67a (5′-AGACAGAAATGGAATAGCACTTTAGACGGGATAAAGGAAAAATTAGCAAGATAC-3′); and for pair 55, M7a (5′-GAACAAAAATGGAGTACAACTTTAAAAAGGGTAGAGAAAAAATTAAAAGAGCAC-3′), M4.1 (5′-GGAAAAGAATGGAATACAACTTTAACAAGGGTAAAGGAAAGATTAAAAAAGCAC-3′), and M2.13 (5′-AGAAAAAAATGGAATATAACTTTAGAAAGGGTAAAGGAAAGATTAAAAGAGCAC-3′).

Open in a separate window

FIG. 3.

Neutralization sensitivities of donor-recipient α2 helix chimeras. (A) An alignment of the predicted amino acid sequences for the α2 helix region is shown for Envs from two transmission pairs (55 [male to female] and 135 [female to male]). The consensus subtype C sequence is shown at the top of each set of sequences. Red arrows show NAb signature residues 335, 336, 337, 343, and 350. The recipient sequences are shown below the consensus (55 F4a, 55 F28a, and 135 M1a). Two recipient Envs with different degrees of neutralization sensitivity were evaluated for pair 55. (B) The neutralization sensitivity of the parental and chimeric Env pseudotypes was evaluated using contemporaneous plasma from donor 55 M or 135 F. The percent virus infectivity relative to the control (lacking donor plasma) is plotted along the vertical axis, and the reciprocal plasma dilution of donor plasma (from 1:20 or 1:100 to 1:62,500) is plotted along the horizontal axis on a log₁₀ scale. Each panel contains infectivity curves for the parental donor (green lines) and recipient (blue lines) Envs and the α2 helix-chimeric recipient Envs (dashed black lines). Infectivity curves for pair 55 and 135 chimeras are shown on the left and right, respectively.

For PCR amplification of the recipient Env (minus the helix) and pCR3.1 vector sequences, 5′-phosphorylated primers were synthesized and gel purified by IDT Technologies. Primer sequences for amplification of the recipient Env backbone (plus vector sequences) and their HXB2 locations were as follows: for pair 135, forward primer 5′-P-CTTCCCTGATAAAATAATAAAATTT-3′ (HXB2 nt 7280 to 7304) and reverse primer 5′-P-ACTAATGTTACAATGTGCTTGTCT-3′ (HXB2 nt 7203 to 7226), and for pair 55, forward primer 5′-P-TTCCCTAATAAAATAATAAAATTT-3′ (HXB2 nt 7281 to 7304) and reverse primer 5′-P-ACTAATGTTACAATGTGCTTGTCT-3′ (HXB2 nt 7203 to 7226).

PCR amplification conditions for the recipient Env backbones were as follows: 1 cycle of 95°C for 3 min; 35 cycles of 95°C for 1 min, 62°C for 30 s, and 72°C for 10 min; 1 cycle of 72°C for 15 min; and storage at 4°C. Pfu Turbo DNA polymerase (Stratagene) was used to generate a blunt-ended PCR amplicon, which was gel purified from an agarose gel by using a QIAquick gel extraction kit (QIAGEN) prior to ligation. The 25-μl PCR mix contained 50 ng of each primer, 10 ng of the plasmid template, 2.5 mM MgCl₂, a 0.2 mM concentration of each deoxynucleoside triphosphate, and 1× reaction buffer. Each synthesized α2 helix DNA fragment was ligated to the purified recipient env backbone by using T4 DNA ligase (5 U/μl; Roche) at 4°C overnight. A portion of the ligation reaction mix was transformed into maximum efficiency XL2-Blue supercompetent cells (1 × 10⁹ CFU/μg DNA; Stratagene). The entire transformation mix was plated onto LB-ampicillin agar plates, generally resulting in ~50 colonies per ligation reaction. These colonies were screened by PCR to identify colonies in which the fragments ligated together in the correct orientation, using forward primer EnvA and a reverse primer that anneals to the 3′ end of the helix. The colony screening primers were as follows: forward, 5′-CAGCACAGTACAATGTACACATGGAA-3′ (HXB2 nt 6950 to 6975); and reverse, 5′-TTCCTTTATCCCGTCTAAAGT-3′ (HXB2 nt 7254 to 7274) for 135 F42a or 135 F67a α2 helix screening, 5′-TGTCTTTACCTCTTCTAAAGT-3′ (HXB2 nt 7254 to 7274) for 135 F48a α2 helix screening, 5′-TTTTAATCTTTCCTTTACCCT-3′ (HXB2 nt 7254 to 7274) for 55 M4.1 or 55 M2.13 α2 helix screening, or 5′-TTTTAATTTTTTCTCTACCCT-3′ (HXB2 nt 7254 to 7274) for 55 M7a α2 helix screening. Colonies that were positive by PCR screening were inoculated into LB-ampicillin overnight cultures, and the plasmids were prepared using a QIAprep Spin miniprep kit (QIAGEN). Plasmids were then screened by (i) restriction digestion, (ii) biological function assay (described below), and (iii) nucleotide sequencing of the α2 helix and flanking regions to confirm that they contained the correct sequences.

Acquisition of length in hypervariable domains is correlated with autologous NAb resistance in donor Envs.

Length polymorphisms in the hypervariable domains V1V2 and V4 were reflected in the donor sequences, where the V1-to-V4 region ranged from 264 to 314 amino acids in length (Fig. 1C and D). The variation in the recipient sequences was less broad, with lengths between 264 and 282 amino acids. Longer hypervariable domains (i.e., increased V1-to-V4 length) were significantly associated with resistance to NAb in linked donor plasma (P = 0.01) (Fig. ​(Fig.1C).1C). This association was recently confirmed for a subset of these Envs, using chimeric recipient Env pseudoviruses in which the native V1V2 domain was replaced with “long” V1V2 domains derived from matched donor Envs (66). In the referenced study, five of seven long donor V1V2 domains tested conferred NAb resistance to the chimeric recipient Env. A significant correlation was not reached between V1-to-V4 length and neutralization by heterologous pooled plasma (Fig. ​(Fig.1D),1D), but there was a trend in this direction (P = 0.08). Thus, hypervariable domain insertions could have the strongest effects on masking strain-specific epitopes, but some protection against cross-neutralization could also be afforded.

Accumulation of N-Gly sites is not correlated with NAb resistance.

Tandem repeat sequences inserted within the longer V1V2 and V4 domains often encoded N-Gly sites (66; data not shown), and in four of the five donors, the number of N-Gly sites was correlated with V1-to-V4 length, using a Spearman rank correlation test (data not shown). Nevertheless, the number of N-Gly sites was not significantly correlated with neutralization by either type of plasma (Fig. 1E and F). This finding could reflect a lack of statistical power for the number of samples analyzed, given the relatively small window of N-Gly sites (17 to 23 sites for all but one Env). Alternatively, the relative positions and packing of the glycans could be more important than net gains or losses for neutralization escape (74).

Individual positions in V1 to V4 are associated with the neutralization phenotype.

A mutual information analysis based on Shannon entropy was applied to the V1-to-V4 sequence and NAb data for the donor and recipient Envs to search for individual positions in the sequence that were predictive of a neutralization-resistant phenotype. Nine “signature” positions were significantly associated with the NAb phenotype after a Bonferroni correction for multiple tests (uncorrected P value, <0.0001) and were mapped three-dimensionally onto a gp120 model that was based on the X-ray structure of the CD4-bound YU2 gp120 core with modeled variable loops (Fig. ​(Fig.2A).2A). Together, the signature residues form a linear display along the immunologically silent face of gp120 (75) and are proximal to the α2 helix (Fig. ​(Fig.2B).2B). Five of the nine residues are located within the α2 helix, and these five positions overlap those previously demonstrated to be under strong positive selection, with dN/dS ratios of >3, using subtype C sequences in the database (Fig. ​(Fig.2B)2B) (22). The convergence of these two independent analyses supports a role for these sequences in neutralization escape. Also, since the database sequences were drawn from a broader sampling that did not include the Envs studied here, these findings could be generally relevant to subtype C viruses and not limited to those circulating in Zambia. Importantly, the positions of the NAb escape signature residues are not subject to the same degree of positive selection in subtype B sequences, where the highest dN/dS ratios are concentrated in V3 (Fig. ​(Fig.2C).2C). Taken together, these data suggest that substitutions in this region could be driven by the autologous NAb response and that changes in this region, in conjunction with alterations in the hypervariable domains, could contribute to neutralization escape during subtype C infection.

Open in a separate window

FIG. 2.

Three-dimensional structural illustration of positions associated with NAb escape and positive selection. Locations were mapped onto a model of gp120 based on the X-ray structure of the CD4-bound YU2 gp120 core (43) (PDB code 1RZK). V1V2 and V3 loops were modeled onto the core for completeness. In the orientation shown, the viral and cellular membranes would be located above and below the protein, respectively. (A) Space-filling model showing locations that were significantly associated with the neutralization phenotype in the mutual information analysis (red balls; HXB2 positions 236, 305, 332, 335, 336, 337, 343, 350, and 393). Positions 335, 336, 337, 343, and 350 are located in the α2 helix. (B and C) Cartoon diagrams showing locations under positive selection, as determined by dN/dS ratios for subtype C sequences (B) and subtype B sequences (C) from the Los Alamos Database (22). Red indicates strong positive selection (dN/dS ratio of >3). Residues associated with the neutralization phenotype in the mutual information analysis are shown as spheres in panel B. In the orientation shown in panels B and C, the viral and cellular membranes would be located above and below the protein, respectively. The orientation in panel A roughly corresponds to rotating the structure in panel B about an axis normal to the viral membrane ~90 degrees counterclockwise.

We gratefully acknowledge the contributions of the staff, lab technicians, participants, and project management group of the Lusaka cohort. We thank Peter Kwong for providing us with the HXB2 trimer coordinates.

This work was supported by NIH grants R01-AI-58706 (C.D.) and R01-AI-51231 (E.H.), by The Bill & Melinda Gates Grand Challenges Program (G.M.S.), and by Yerkes National Primate Research Center base grant RR-00165, awarded by the National Center for Research Resources of the National Institutes of Health.