1. Evolutionary Biology
  2. Microbiology and Infectious Disease
Download icon

Mapping mutational effects along the evolutionary landscape of HIV envelope

  1. Hugh K Haddox
  2. Adam S Dingens
  3. Sarah K Hilton
  4. Julie Overbaugh
  5. Jesse D Bloom  Is a corresponding author
  1. Fred Hutchinson Cancer Research Center, United States
  2. University of Washington, United States
Research Article
  • Cited 4
  • Views 1,600
  • Annotations
Cite this article as: eLife 2018;7:e34420 doi: 10.7554/eLife.34420

Abstract

The immediate evolutionary space accessible to HIV is largely determined by how single amino acid mutations affect fitness. These mutational effects can shift as the virus evolves. However, the prevalence of such shifts in mutational effects remains unclear. Here, we quantify the effects on viral growth of all amino acid mutations to two HIV envelope (Env) proteins that differ at >100 residues. Most mutations similarly affect both Envs, but the amino acid preferences of a minority of sites have clearly shifted. These shifted sites usually prefer a specific amino acid in one Env, but tolerate many amino acids in the other. Surprisingly, shifts are only slightly enriched at sites that have substituted between the Envs—and many occur at residues that do not even contact substitutions. Therefore, long-range epistasis can unpredictably shift Env’s mutational tolerance during HIV evolution, although the amino acid preferences of most sites are conserved between moderately diverged viral strains.

https://doi.org/10.7554/eLife.34420.001

eLife digest

The virus that causes AIDS, or HIV, has a protein called Env on its surface, which is essential for the virus to infect cells. Env can also be recognized by the immune system, which then targets the virus for destruction or blocks it from infecting cells. Unfortunately, Env evolves very quickly, which means that HIV can evade our defenses. However, there are limits to how much this protein can change, since it still needs to perform its essential role in helping viruses enter cells.

In the century since HIV first appeared in human populations, the virus has evolved considerably. There are now many HIV strains that infect people, and they bear Env proteins with substantially different sequences. However, it is not clear if these changes in sequence have resulted in Envs from distinct strains being able to tolerate different mutations.

To examine this question, Haddox et al. compared how the Envs from two strains of HIV react to modifications in their sequences. They created all possible individual mutations in the proteins, and the resulting collections of mutated viruses were then tested for their ability to infect cells in the laboratory.

Most mutations had similar effects in both Env proteins. This allowed Haddox et al. to identify portions of the protein that easily accommodate changes, and portions that must remain unchanged for viruses to remain infectious—at least in the laboratory. Some of these mutations are under different types of pressures when the virus faces the immune system, and those were identified using computational approaches.

However, some mutations were tolerated differently by the two Env proteins. Therefore, viral strains differ in how their Env proteins can evolve. The parts of Env that showed differences in mutational tolerance between the strains were not necessarily the parts that differ in sequence. This shows that changes in sequence in one part of the protein can modify how other portions evolve.

It remains to be determined whether changes in tolerance to mutations translate into differences in how the virus can escape immunity. This is an important question given that the rapid evolution of Env is a major obstacle to creating a vaccine for HIV.

https://doi.org/10.7554/eLife.34420.002

Introduction

HIV’s envelope (Env) protein evolves very rapidly. The major group of HIV-1 that is responsible for the current pandemic originated from a virus that entered the human population 100 years ago (Sharp and Hahn, 2011; Worobey et al., 2008; Faria et al., 2014). The descendants of this virus have evolved so rapidly that their Envs now have as little as 65% protein identity (Lynch et al., 2009). For comparison, protein orthologs shared between humans and mice have only diverged to a median identity of 78% over 90 million years (Waterston et al., 2002; Hedges et al., 2006).

Env’s rapid evolution has dire consequences for anti-HIV immunity, since it erodes the efficacy of most neutralizing antibodies (Albert et al., 1990; Wei et al., 2003; Richman et al., 2003; Burton et al., 2005). Because of this public-health importance, numerous studies have experimentally characterized aspects of the ‘evolutionary landscape’ that Env traverses. The immediate evolutionary space accessible to any given Env is largely defined by the effects on viral fitness of all single amino acid mutations to Env. Most mutational studies have measured how just a small number of these mutations affect viral growth in cell culture, although it has recently become possible to use deep mutational scanning to measure the effects of many (Al-Mawsawi et al., 2014; Duenas-Decamp et al., 2016) or even all (Haddox et al., 2016) single amino acid mutations to an Env variant.

But interpreting these studies in the context of Env evolution requires addressing a fundamental question: How informative are mutational studies of a single protein variant about constraints on long-term evolution? During protein evolution, substitutions at one site can change the effect of mutations at other sites (Natarajan et al., 2013; Gong et al., 2013; Harms and Thornton, 2014; Podgornaia and Laub, 2015; Starr and Thornton, 2016; Klink and Bazykin, 2017). We will follow the nomenclature of (Pollock et al., 2012) to refer to these changes in mutational effects as shifts in a site’s amino acid preferences. Such shifts can accumulate as substitutions become entrenched via epistatic interactions with subsequent changes (Starr et al., 2017; Pollock et al., 2012; Shah et al., 2015; Bazykin, 2015)—although the magnitude of these shifts is usually limited (Doud et al., 2015; Chan et al., 2017; Ashenberg et al., 2013; Risso et al., 2015).

Given that the Envs of circulating HIV strains represent a vast collection of homologs that often differ at >100 residues, shifts in amino acid preferences could make the outcome of any study highly dependent on the Env used. Indeed, epistasis among a few combinations of Env mutations has been experimentally demonstrated (da Silva et al., 2010), and epistatic fitness landscapes have been computationally inferred for a variety of HIV proteins (Kouyos et al., 2012; Ferguson et al., 2013; Mann et al., 2014; Barton et al., 2015) including Env (Louie et al., 2018). However, the only protein-wide experimental studies of how amino acid preferences shift during evolution have examined proteins that are structurally far simpler than Env, which forms a large heavily glycosylated heterotrimeric complex that transitions through multiple conformational states (Munro et al., 2014; Ozorowski et al., 2017).

Here, we use an improved version of a previously described deep mutational scanning strategy (Haddox et al., 2016) to measure the effects on viral growth of all single amino acid mutations to two transmitted-founder virus Envs that differ by >100 mutations. We compare these complete maps of mutational effects to identify sites that have shifted in their amino acid preferences between the Envs. Most sites show no detectable shifts, but 30 sites have clearly shifted preferences. These shifted sites usually prefer a specific amino acid in one Env but have shifted to tolerate many amino acids in the other Env. The shifted sites cluster in structure but are often distant from any amino acid substitutions that distinguish the two Envs, demonstrating the action of long-range epistasis. By aggregating our measurements for both Envs, we identify sites that evolve faster or slower in nature than expected given the functional constraints measured in the lab, probably due to pressure for immune evasion. Overall, our work provides complete across-strain maps of mutational effects that inform analyses of Env’s evolution and function.

Results

Two Envs from clade A transmitted-founder viruses

The viruses most relevant to HIV’s long-term evolution are those which are transmitted from human-to-human. However, the only prior work that has measured how all Env amino acid mutations affect HIV growth is a study by some of us (Haddox et al., 2016) that used a late-stage lab-passaged CXCR4-tropic virus (Peden et al., 1991). The properties of Env can vary substantially between such late-stage viruses and the transmitted-founder viruses relevant to HIV’s long-term evolution (Sagar et al., 2006; Wilen et al., 2011; Parrish et al., 2013; Ronen et al., 2015).

We therefore selected Envs from two transmitted-founder viruses, BG505.W6M.C2.T332N and BF520.W14M.C2 (hereafter referred to as BG505 and BF520), that were isolated from HIV-infected infants shortly after mother-to-child transmission (Nduati et al., 2000; Wu et al., 2006; Goo et al., 2014). The BG505 Env has been extensively studied from a structural standpoint (Julien et al., 2013; Lyumkis et al., 2013; Pancera et al., 2014; Huang et al., 2014; Sanders et al., 2015; Stewart-Jones et al., 2016; Gristick et al., 2016), and variants of this Env are being tested as vaccine immunogens (Sanders et al., 2013, 2015; de Taeye et al., 2015). We used the T332N variant of BG505 Env because it has a common glycosylation site that is targeted by many anti-HIV antibodies (Sanders et al., 2013). The BF520 Env was isolated from an infant who developed an early broad anti-HIV antibody response (Goo et al., 2014; Simonich et al., 2016). We have previously created comprehensive codon-mutant libraries of the BF520 Env and used them to map HIV antibody escape (Dingens et al., 2017), but these BF520 libraries have not been characterized with respect to how mutations affect viral growth.

Both BG505 and BF520 are from clade A of the major (M) group of HIV-1. Figure 1 shows the phylogenetic relationship among these two Envs and other clade A sequences. BG505 and BF520 are identical at 721 of the 836 pairwise-alignable protein sites (86.2% identity). However, in our experiments we mutagenized only the ectodomain and transmembrane domain of Env, and excluded the signal peptide and cytoplasmic tail. The reason is that we measure how Env mutations affect viral growth, which is influenced both by the functionality of Env protein molecules and their expression level. Mutations in the signal peptide and cytoplasmic tail commonly affect Env expression level (Chakrabarti et al., 1989; Yuste et al., 2004; Li et al., 1994), so we excluded these regions with the goal of reducing the degree to which we simply identified mutations that affected Env expression. In the ectodomain and transmembrane domains of Env, BG505 and BF520 are identical at 549 of the 616 sites (89.1% identity) that are alignable across clade A Envs (Figure 1—source data 1, Figure 1—source data 2). The divergence between BG505 and BF520 therefore offers ample opportunity to investigate mutational shifts during Env evolution.

Figure 1 with 1 supplement see all
Phylogenetic tree showing the relationship of BG505 and BF520 to other clade A Envs.

The tree shows the 69 Envs in the alignment in Figure 1—source data 1, which is a subsample of clade A sequences from the group M alignment in the Los Alamos HIV sequence database (http://www.hiv.lanl.gov). Sites not mutagenized in our experiments (the signal peptide and cytoplasmic tail) or that are poorly alignable were masked as indicated in Figure 1—source data 2, leaving 616 alignable sites. The pairwise identity of BG505 and BF520 to other sequences at alignable sites is in Figure 1—figure supplement 1. The tree topology was inferred using RAxML (Stamatakis, 2014) under the GTRCAT model of nucleotide substitution, and branch lengths were optimized under the M0 Goldman-Yang model (Yang et al., 2000) using phydms (Hilton et al., 2017).

https://doi.org/10.7554/eLife.34420.003

Deep mutational scanning of each Env

We have previously described a deep mutational scanning strategy for measuring how all amino acid mutations to Env affect HIV growth in cell culture, and applied this strategy to the late-stage lab-adapted LAI strain (Haddox et al., 2016). Here, we made several modifications to this earlier strategy to apply it to transmitted-founder Envs and to reduce the experimental noise. This last consideration is especially important when comparing Envs, since it is only possible to reliably detect differences that exceed the magnitude of the experimental noise. Our modified deep mutational scanning strategy is in Figure 2A. This approach had the following substantive changes: instead of SupT1 cells, we used SupT1.CCR5 cells (SupT1 cells that express CCR5 in addition to CXCR4 [Boyd et al., 2015]) to support growth of viruses with transmitted-founder, CCR5-tropic Envs; we used more virions for the first passage (3×106 versus 5×105 infectious units per library) to avoid bottlenecking library diversity; and rather than performing a full second passage we just did a short high-MOI infection to enable recovery of env genes from infectious virions without bottlenecking (Figure 2A). We performed this deep mutational scanning in full biological triplicate for both BG505 and BF520 (Figure 2B). Our libraries encompassed all codon mutations to all sites in Env except for the signal peptide and cytoplasmic tail.

Figure 2 with 1 supplement see all
Deep mutational scanning workflow.

(A) We made libraries of proviral HIV plasmids with random codon-level mutations in the env gene. The number of mutations per gene approximately followed a Poisson distribution with a mean between 1 and 1.5 (Figure 2—figure supplement 1). We transfected the plasmids into 293T cells to generate mutant viruses, which lack a genotype-phenotype link since cells are multiply transfected. To establish a genotype-phenotype link and select for Env variants that support HIV growth, we passaged the libraries in SupT1.CCR5 cells for four days at a low multiplicity of infection (MOI) of 0.01. To isolate the env genes from only viruses that encoded a functional Env protein, we infected the passaged libraries into SupT1.CCR5 cells at high MOI and harvested reverse-transcribed non-integrated viral DNA after 12 hr. We then deep sequenced the env genes from these final samples as well as the initial plasmid library, using molecular barcoding to reduce sequencing errors. We also deep sequenced identically handled wildtype controls to estimate error rates. Using these sequencing data, we estimated the preference for each of the 20 amino acids at each site in Env. These data are represented in logo plots, with the height of each letter proportional to that site’s preference for that amino acid. (B) We conducted this experiment in full biological triplicate for both BG505 and BF520, beginning each replicate with independent creation of the plasmid mutant library. These replicates therefore account for all sources of noise and error in the experiments.

https://doi.org/10.7554/eLife.34420.007

The deep mutational scanning effectively selected for functional Envs as evidenced by strong purifying selection against stop codons. Figure 3A shows the average frequency of mutations across Env in the plasmid mutant libraries, the mutant viruses, and wildtype controls as determined from the deep sequencing. The mutant viruses show clear selection against stop codons and many nonsynonymous mutations (Figure 3A). This selection is more apparent if we correct for the background error rates estimated from the wild-type controls (Figure 3—source data 1). The error-corrected frequencies of stop codons drop to 3–16% of their original values (Figure 3—source data 1), with the residual stop codons probably due to some non-functional virions surviving due to complementation by other co-infecting virions. The error-corrected frequencies of nonsynonymous mutations also drop substantially (43%–49% of their original values), whereas the frequencies of synonymous mutations drop only slightly (85%–95% of their original values). These trends are consistent with the fact that nonsynonymous mutations are often deleterious, whereas synonymous mutations often (Zanini and Neher, 2013) have only mild effects on viral growth. Figure 3A only summarizes one aspect of the deep mutational scanning data, but Supplementary file 1 and 2 contain detailed plots showing all aspects of the data (read depth, per-site mutation rate, etc) as generated by the dms_tools2 software (Bloom, 2015https://jbloomlab.github.io/dms_tools2/).

Figure 3 with 1 supplement see all
The deep mutational scanning selects for functional Envs and yields measurements that are well correlated among replicates.

(A) The average per-codon mutation frequency when sequencing plasmids encoding wildtype Env (DNA), plasmid mutant libraries (mutDNA), mutant viruses after the final infection (mutvirus), and virus generated from wild-type plasmids (virus). Mutations are categorized as nonsynonymous, synonymous, or stop codon. The DNA samples show that sequencing errors are rare, and the virus samples show that viral-replication errors are well below the frequency of mutations in the mutDNA samples. Comparing the mutvirus to mutDNA shows clear purifying selection against stop codons and some nonsynonymous mutations, particularly after subtracting the background error rates given by the virus and DNA samples (Figure 3—source data 2). More extensive plots from the analysis of the deep sequencing data are in Supplementary file 1 and 2. (B) Correlations between replicates in the measured preferences of each site in Env for all 20 amino acids. Blue indicates replicate measurements on BF520, red indicates replicate measurements on BG505, and gray indicates across-Env measurements of BF520 versus BG505. R is the Pearson correlation coefficient. The numerical values for the preferences are in Figure 3—source data 2. Figure 3—figure supplement 1 shows the correlations using contour rather than scatter plots.

https://doi.org/10.7554/eLife.34420.009

We used the deep mutational scanning data to estimate the preference of each site in Env for each amino acid via the analysis method described in Bloom (2015). As graphically illustrated in Figure 2A, the preferences for each site are normalized to sum to one. Our libraries were mutagenized at 670 sites in BG505 and 662 sites in BF520, so 670×20=13,400 and 662×20=13,240 preferences were estimated for each Env, respectively. The correlations between the preferences from different experimental replicates are in Figure 3B, and the preferences themselves are in Figure 3—source data 2. These replicate-to-replicate correlations are substantially higher than those for the deep mutational scanning of LAI Env by (Haddox et al., 2016), which had replicate-to-replicate Pearson correlations of only R=0.45 to 0.50.

While the replicates are well correlated across all replicates for both BG505 and BF520, the replicates for BG505 are more correlated with each other than with replicates for BF520, and vice versa (Figure 3B, compare red and blue versus gray plots). This fact hints that there are some shifts in amino-acid preferences between the two Envs—something that is investigated with more statistical rigor later in this paper. Note also that there is a trend for highly preferred amino acids to be more strongly preferred in BG505 than BF520 (most high-preference points in the gray plots in Figure 3B fall above the diagonal); however, this trend does not necessarily reflect differences between the Envs. Rather, there were modest differences in the stringency of selection between our deep mutational scans of BG505 and BF520 (Figure 3—source data 1 shows that purifying selection better purged stop codons in BG505). In the next section, we correct for these experimental differences by calibrating each dataset to match the stringency of selection in nature.

Amino acid preferences of the Envs and their relationship to HIV evolution

The most immediate question is how authentically the experimental measurements describe the actual selection on Env function in nature. Direct comparisons between experimentally measured amino acid preferences and amino acid frequencies in natural sequences are confounded by the fact that the natural sequences are evolutionarily related. This problem can be overcome by making the comparison in a phylogenetic context to account for the evolutionary relationships among sequences.

Specifically, we used our deep mutational scanning data to construct experimentally informed codon models (ExpCM’s) for Env’s evolution. An ExpCM is a phylogenetic substitution model that incorporates the functional constraints measured in a deep mutational scanning experiment (Hilton et al., 2017). If the experiment captures much of the actual evolutionary constraint on a gene, then an ExpCM will describe the gene’s natural evolution better than a standard phylogenetic codon substitution model. The reason is that standard codon substitution models (Yang et al., 2000) only model functional constraint via a single parameter that represents the rate of fixation of nonsynonymous protein-altering mutations relative to synonymous ones; this parameter is called dN/dS or ω. In contrast, an ExpCM accounts for the preference of each site for each of the 20 amino acids under the functional selection in the deep mutational scan, and then additionally adds an ω parameter that represents the relative rate of nonsynonymous to synonymous substitutions after accounting for these functional constraints (Bloom, 2017; Hilton et al., 2017). Importantly, since we expect some sites in Env to be under diversifying selection from immunity, we extended the ExpCM’s described in Hilton et al. (2017) to draw ω from a gamma distribution as is commonly done for codon-substitution models (Yang et al., 2000).

Table 1 shows that ExpCM’s informed by the deep mutational scanning of either BG505 or BF520 describe the natural evolution of Env vastly better than a standard codon substitution model. In addition to the improved fit of the ExpCM’s, we can also interpret the ω parameter. Recall that for standard codon substitution models, ω is simply the rate of fixation of nonsynonymous mutations relative to synonymous ones. For such models, the gene-wide average ω is almost always <1, since purifying selection purges many functionally deleterious amino acid mutations even for adaptively evolving proteins (Murrell et al., 2015). Indeed, Table 1 shows that Env’s gene-wide average ω is <one for a standard model. But for ExpCM’s, ω is the relative rate of nonsynonymous to synonymous substitutions after accounting for functional constraints measured in the deep mutational scanning (Bloom, 2017). For the ExpCM’s, the gene-wide average ω is >1 (Table 1), indicating that external selection (e.g. from immunity) drives Env to fix amino acid mutations faster than expected under a null model that only accounts for functional constraints on the protein.

Table 1
Evolutionary models informed by the deep mutational scanning describe HIV’s evolution in nature much better than a standard substitution model.

Shown are the results of maximum likelihood fitting of substitution models to the clade A phylogeny in tree. Experimentally informed codon models (Hilton et al., 2017) utilizing the across-replicate average of the deep mutational scanning describe Env’s natural evolution far better than a standard codon substitution model (Yang et al., 2000) as judged by comparing the Akaike information criteria (Posada and Buckley, 2004). Both ExpCM’s have a stringency parameter >1. All models draw ω from a gamma distribution, and the table shows the mean (ω¯) and shape parameters (ωα and ωβ) of this distribution. The last two columns show the number of sites evolving faster (ωr>1) or slower (ωr<1) than expected at a false discovery rate of 0.05, as determined using the approach in Bloom (2017) (see also the last section of the Results). Analyses were performed using phydms (Hilton et al., 2017, http://jbloomlab.github.io/phydms/). Table 1—source data 1. shows the results for additional substitution models.

https://doi.org/10.7554/eLife.34420.013
ModelΔAICLogLikelihoodnParamsStringency

ω¯

ωα

ωβ

Nsites ωr>1Nsites ωr<1
ExpCM BF5200.0−35218.872.81.41.00.76635
ExpCM BG505269.0−35353.372.11.30.90.76553
Goldman-Yang M53455.1−36941.412nan0.80.60.714211

ExpCM’s also have a stringency parameter that relates selection in the experiments to that in nature. Essentially, this parameter indicates how strongly natural selection prefers the amino acids that are preferred in the deep mutational scanning (Hilton et al., 2017). A stringency parameter >1 indicates that natural selection prefers the same amino acids as the experiments, but with greater stringency. Both ExpCM’s have stringency parameters >1 (Table 1)—a finding that makes sense, since the stop-codon analysis in the previous section suggests that the experimental selections are more lax than natural selection on HIV.

For the entire rest of the paper, we use the experimentally measured preferences re-scaled by the stringency parameters in Table 1. The reason we do this is to distinguish genuine differences between the two Envs from mere variation in the strength of selection between the two sets of experiments. Re-scaling both sets of preferences to optimally describe Env evolution in nature is a principled way to standardize the measurements; see (Hilton et al., 2017) and the Materials and methods section entitled ‘Re-scaling the preferences’ for a more detailed explanation.

A qualitative way to assess if the deep mutational scanning authentically describes selection on Env function is to visually compare the measurements with existing knowledge. Figure 4 and Figure 5 show the re-scaled across-replicate average of the amino acid preferences for each Env. At sites of known functional importance, these preferences are usually consistent with prior knowledge. For instance, residues T257, D368, E370, W427, and D457 are important for Env binding to CD4 (Olshevsky et al., 1990), and all these amino acids are highly preferred in our deep mutational scanning (Figure 4 and Figure 5). Likewise, Env has 10 disulfide bonds (linking sites 54–74, 119–205, 126–196, 131–157, 218–247, 228–239, 296–331, 378–445, 385–418, and 598–604), most of which are important for function (van Anken et al., 2008)—and the cysteines at these sites are highly preferred in our deep mutational scanning. The deep mutational scanning is also consistent with prior knowledge about sites that are tolerant of mutations. For instance, Env has five variable loops that mostly evolve under weak constraint in nature (Starcich et al., 1986; Zolla-Pazner and Cardozo, 2010)—and most sites in these loops are mutationally tolerant in our deep mutational scanning (see sites indicated by gray overlay bars in Figure 4 and Figure 5, such as 132 to 195). It is beyond the scope of this paper to catalog associations between our measurements and all other prior mutational studies of Env, but the concordance of our findings with the above mutational studies, and the fact that our data improve phylogenetic models of Env’s natural evolution, suggest that our experiments do a reasonable job of authentically measuring functional selection on Env.

Amino acid preferences for the BG505 Env.

At each site, the height of the letter is proportional for that site’s preference for that amino acid. The top color bar indicates the region of Env (gp120 variable loop, gp120 not variable loop, or gp41). The lower color bar indicates the evidence that the site evolves faster (ωr>1) or slower (ωr<1) than expected given the experiments (Bloom, 2017). We report the p-value for ωr1 rather than the value of ωr itself since point estimates of ωr are unreliable for individual sites due to low numbers of observations, making the p-value a better indicator of the strength of the statistical evidence for faster or slower than expected evolution (Kosakovsky Pond and Frost, 2005; Murrell et al., 2012). The letters above the logos indicate the wildtype amino acid in BG505. Sites are numbered using the HXB2 scheme (Korber et al., 1998). This logo plot shows the site-specific amino acid preferences for BG505 after averaging the replicates and re-scaling by the stringency parameter in Table 1. The figure was generated using dms_tools2 (Bloom, 2015), which in turn utilizes weblogo (Crooks et al., 2004). The numerical values of the preferences are in Figure 4—source data 1, the mapping from sequential to HXB2 numbering is in Figure 4—source data 2, and the ωr values are in Figure 4—source data 3.

https://doi.org/10.7554/eLife.34420.015
Amino acid preferences for the BF520 Env.

This figure is the same as Figure 4 except that it shows the data for BF520 instead of BG505. The numerical values of the preferences are in Figure 5—source data 1, the mapping from sequential to HXB2 numbering is in Figure 5—source data 2, and the ωr values are in Figure 5—source data 3.

https://doi.org/10.7554/eLife.34420.019

Shifts in amino acid preferences between BG505 and BF520

The most fundamental question that we seek to address is how similar the amino acid preferences are between the two Envs. We have already noted that Figure 3B shows that the preferences are more correlated for replicate measurements on the same Env than for replicate measurements on different Envs. However, simply comparing correlation coefficients does not identify specific sites where mutational effects have shifted, nor does it quantify the magnitude of any shifts.

We therefore used a more rigorous approach to identify sites where the amino acid preferences differ between BG505 and BF520 by an amount that exceeds the noise in our experiments. We first re-scaled the preferences from each experimental replicate by the stringency parameter for that Env from Table 1 to calibrate all measurements to the stringency of natural selection. We then identified the 659 sites in the mutagenized regions of Env that are pairwise alignable between BG505 and BF520 (Figure 6—source data 1). For each site, we calculated the shift in amino acid preferences between Envs using an approach similar to that of (Doud et al., 2015) as illustrated in Figure 6A. This approach calculates the magnitude of the shift after correcting for experimental noise by comparing the differences in preferences between replicates for BG505 and BF520 to the differences between replicates for the same Env. Figure 6A shows this calculation for a site that has not shifted (site 598, which strongly prefers cysteine in both Envs), the most shifted site (512, which shifts from being mutationally tolerant in BG505 to strongly preferring alanine in BF520), and two other sites with more intermediate behaviors.

Env sites with shifted amino acid preferences between BG505 and BF520.

Note that the preferences have been re-scaled using the stringency parameters in Table 1 to enable direct comparison across Envs. (A) Calculation of the corrected distance between the amino acid preferences of BG505 and BF520 at four example sites. We have triplicate measurements for each Env. We calculate the distance between each pair of replicate measurements, and group these into comparisons between the two Envs and within replicates for the same Env. We compute the root-mean-square distance (RMSD) for both sets of comparisons, which we denote as RMSDbetween and RMSDwithin. The latter quantity is a measure of experimental noise. The noise-corrected distance between Envs at a site, RMSDcorrected, is simply the distance between the two Envs minus this noise. (B) The bottom distribution (orange) shows the corrected distances between BG505 and BF520 at all alignable sites (see Figure 6— source data 1 for numerical values). The next distribution (blue) is a null generated by computing the corrected distances on all randomizations of the replicates among Envs. The top two distributions (green) compare Env to the non-homologous influenza hemagglutinin (HA) protein (Doud and Bloom, 2016) simply putting sites into correspondence based on sequence number. We compute the p-value that a site has shifted between BG505 and BF520 as the fraction of the null distribution that exceeds that shift, and identify significant shifts at a false discovery rate (FDR) of 0.1 using the method of (Benjamini and Hochberg, 1995). Using this approach, 30 of the 659 sites have significant shifts (corrected distance 0.22). (C) All sites that have significantly shifted their amino acid preferences at an FDR of 0.01. For each site, the logo stacks show the across-replicate average preferences for BG505 and BF520. The wild-type amino acid for that Env is indicated using the small black letters above each logo plot; note how the wild-type amino acid is frequently but not always the most preferred one. The sites are sorted by the magnitude of the shift.

https://doi.org/10.7554/eLife.34420.023

The overall distribution of shifts between BG505 and BF520 is shown in Figure 6B. Most sites have relatively small shifts (close to zero), although there is a long tail of sites with large shifts. This tail reaches its upper value with site 512, which has a shift of 0.52 out of a maximum possible of 1.0. How should we interpret this distribution—have mutational effects shifted a lot, or not very much? We can establish an upper-bound for how much sites might shift by comparing Env to a non-homologous protein. Figure 6B shows the distribution of shifts when comparing Env to influenza’s hemagglutinin protein, which has previously had its amino acid preferences measured by deep mutational scanning (Doud and Bloom, 2016). Most sites have large shifts between Env and hemagglutinin, with the typical shift being 0.4 and some approaching the maximum value of 1.0. We can also establish a lower-bound by creating a null distribution for the expected shifts if all differences are simply due to experimental noise. This null distribution is created by randomizing the experimental replicates among Envs. Figure 6B shows that the null distribution is more peaked at zero than the real distribution, and does not have the same prominent tail of sites with large shifts. The answer to the question of how much mutational effects have shifted is therefore nuanced: they have substantially shifted at some sites, but remain vastly more similar between the two Envs than between two unrelated proteins.

We can use the null distribution to identify sites where the shifts between BG505 and BF520 are significantly larger than the noise in our experiments (Figure 6B). There are 30 such sites at a false discovery rate of 0.1. Figure 6C shows the amino acid preferences of these significantly shifted sites for each Env. For the majority of shifted sites, one Env prefers a specific amino acid whereas the other Env tolerates many amino acids; for instance, see sites 512, 516, 599, 165, 605 and 505 in Figure 6C. Such broadening and narrowing of a site’s mutational tolerance is frequently linked to changes in protein stability, with a more stable protein typically being more mutationally tolerant (Wang et al., 2002; Bloom et al., 2006; Gong et al., 2013; Kumar et al., 2017). Work with engineered Env protein in the form of ‘SOSIP’ trimer (Binley et al., 2000; Sanders et al., 2002) has shown that BG505 SOSIP is more thermostable than BF520 SOSIP (Verkerke et al., 2016). Consistent with this fact, sites with altered mutational tolerance are often (although not always, see sites 165 and 520 in Figure 6C) more mutationally tolerant in BG505. Differences in Env’s expression level might also contribute to a general broadening or narrowing of tolerance to subsequent mutations. The reason is that our experiments select for viral growth (which is affected by both Env function and expression), so it is possible that some of the shifts are due to epistatic mutational effects on expression rather than function.

However, not all of the significantly shifted sites show a simple pattern of broadening or narrowing of mutational tolerance. For instance, site 288 does not alter its mutational tolerance but rather flips its rather narrow amino acid preference from phenylalanine in BG505 to leucine in BF520 (Figure 6C). Thus, there is variation in both the extent and types of shifts observed.

Structural and evolutionary properties of shifted sites

What distinguishes the sites that have undergone significant shifts? First, we analyzed the distribution of shifted sites in context of Env’s three-dimensional structure. Env’s structure is highly conformationally dynamic and undergoes large changes upon receptor binding and membrane fusion. In an effort to account for these dynamics, we examined multiple conformational states of Env: the closed pre-fusion state (Stewart-Jones et al., 2016), the open CD4-bound state (Ozorowski et al., 2017), and the post-fusion six-helix bundle (Weissenhorn et al., 1997). Figure 7A shows the locations of the shifted sites on the crystal structure of Env in the closed pre-fusion state. There is no visually obvious tendency for shifted sites to preferentially be on Env’s surface or in its core, and statistical analysis of both the closed and open states of Env (Figure 7B) finds no association between a site’s relative solvent accessibility and whether its amino acid preferences have shifted. We did not attempt to analyze the association between solvent accessibility and shift for the post-fusion six-helix bundle because crystal structures of this conformation only contain 80 Env residues (Weissenhorn et al., 1997; Chan et al., 1997; Tan et al., 1997).

Figure 7 with 2 supplements see all
Characteristics of significantly shifted sites.

(A) One monomer of the closed pre-fusion Env trimer (Stewart-Jones et al., 2016) is colored from white to red according to the magnitude of the mutational shift at each site (red indicates large shift). Sites that are significantly shifted according to Figure 6B are in spheres, and all other sites are in cartoon representation. (B) There is no significant difference in the relative solvent accessibility of sites that have and have not undergone significant shifts. This observation holds for both the closed trimer conformation in (A) and the CD4-bound trimer conformation (Ozorowski et al., 2017). The absolute solvent accessibility of each site was calculated using DSSP (Kabsch and Sander, 1983) and normalized to a relative solvent accessibility using the absolute accessibilities from Tien et al. (2013). (C) Sites of significant shifts are clustered in the structures of both the closed and open Env trimers. The left two plots show the distance of each significantly shifted and not-shifted site to the closest other shifted site in the indicated structure. The right-most plot shows the minimum distance across both conformations. The trend for shifts to cluster becomes stronger when considering the minimum distance, suggesting multiple conformations contribute to this trend. (D) Large mutational shifts are not strongly enriched at sites that have substituted between BG505 and BF520, or at sites that contact sites that have substituted. The plots show the magnitudes of the shifts among structurally resolved sites that have substituted between BG505 and BF520, the non-substituted sites that physically contact a substitution in the indicated structure(s) (any non-hydrogen atom within 3.5 angstroms), and all other sites. Figure 7— source data 1 shows that there is a borderline-significant tendency of significantly shifted sites to have substituted. All plots only show sites that are structurally resolved in the indicated structure(s). Structural distances and solvent accessibilities were calculated using all monomers in the trimer. P-values were calculated using the Mann-Whitney U test. Figure 7—figure supplement 1 and Figure 7—figure supplement 2 zoom in on some relevant clusters of sites.

https://doi.org/10.7554/eLife.34420.025

However, Figure 7A does suggest that the sites of significant shifts tend to cluster in Env’s structure. A statistical analysis confirms that there is clustering of shifted sites for the closed and open conformations, with the effect being strongest when we define contacts based on the closest intra-residue distance across these two conformations (Figure 7C). Therefore, the factors that drive shifts in Env’s mutational tolerance often affect physically interacting clusters of residues in a coordinated fashion. We also investigated clustering of shifted sites in the post-fusion six-helix bundle. Because structures of this conformation only resolve the coordinates of 80 residues, we did not perform a statistical analysis. However, a qualitative analysis revealed that three of the four shifted sites that are resolved in the post-fusion conformation cluster at one end of the helical bundle (Figure 7—figure supplement 1).

An obvious hypothesis is that strongly shifted sites have substituted between BG505 and BF520, or physically contact such substitutions. According to this hypothesis, substitutions would alter the local physicochemical environment of the substituted site and its neighbors, thereby shifting the amino acid preferences of sites in the physical cluster. But surprisingly, for both the closed and open conformations, the typical magnitude of shifts is not significantly larger at sites that have substituted, or at sites that contact sites that have experienced substitutions (Figure 7C). For the six-helix bundle, there are five structurally resolved substituted sites, one of which is adjacent to the cluster of shifted sites (Figure 7—figure supplement 1). The number of resolved shifted and substituted sites in this structure is too small for a meaningful statistical analysis of the type in Figure 7D. However, the cluster of shifted and substituted sites in the six-helix bundle is also present in the closed and open states (Figure 7—figure supplement 1), and so is included in the statistical analyses in Figure 7D.

There is a borderline trend for the significantly shifted sites to be more likely to have substituted between BG505 and BF520 (Figure 7—source data 1), but most shifted sites have not substituted (only 8 of the 30 shifted sites differ in amino acid identity between the two Envs). The lack of strong enrichment in shifts at substituted sites contrasts with previous protein-wide experimental (Doud et al., 2015) and simulation-based (Pollock et al., 2012; Shah et al., 2015) studies of shifting amino acid preferences, which found that shifts were dramatically more pronounced at sites of substitutions. The difference may arise because these earlier studies examined proteins that are fairly conformationally static (absolutely so in the case of the simulations). The fact that Env is extremely complex and conformationally dynamic (Munro et al., 2014; Ozorowski et al., 2017) may increase the opportunities for long-range epistasis to enable substitutions at one site to shift the amino acid preferences of distant sites.

Indeed, many of the shifted sites cluster within regions of Env that are highly conformationally dynamic. Figure 7—figure supplement 2 shows the structural context of these clusters in finer detail. One cluster is at the trimer apex where two of Env’s variable loops pack against one another and against an adjacent protomer. These interactions are likely involved in regulating the transition between conformational states, and upon CD4 binding, these loops become highly disordered (Guttman et al., 2014; Ozorowski et al., 2017). Mutations at two of the shifted sites in this cluster (165 and 307) have been shown to cause Env to assume aberrant conformations, suggesting that these sites can strongly modulate Env’s dynamics (Lee et al., 2017). Strikingly, this cluster of shifted sites may reflect previously observed differences in the conformational dynamics of this regions between these two Envs; the V2 region of BF520 SOSIP trimer is more accessible to deuterium exchange than the BG505 SOSIP trimer (Verkerke et al., 2016). The other cluster of shifted sites is near a network of hydrophobic amino acids that has been proposed to help transmit the large-scale conformational change that takes place upon CD4 binding (Ozorowski et al., 2017). One of the shifted sites (site 69) overlaps with this network, and mutations at another (site 64) have been shown to strongly modulate the relative stability of the open and closed conformations (de Taeye et al., 2015). In total, these two clusters consist of nearly half of the shifted sites (13 out of 30). One hypothesis why so many shifted sites cluster in these regions is that their dynamic nature allows long-range epistatic interactions to be readily propagated between substituted sites and distant shifted sites. It is difficult to discern exactly how these interactions might occur, but there is certainly a trend for sites that are conformationally dynamic to also be sites that show shifts in their amino acid preferences during evolution.

Entrenchment of substitutions modestly contributes to mutational shifts

One idea that has recently gained support in the protein-evolution field is that substitutions become ‘entrenched’ by subsequent evolution (Pollock et al., 2012; Shah et al., 2015; Starr et al., 2017). Entrenchment is the tendency of a mutational reversion to become increasingly unfavorable as a sequence evolves. Given two homologs, if there is no entrenchment then the effect of mutating a site in the first homolog to its identity in the second will simply be the opposite of mutating the site in the second homolog to its identity in the first. But if there is entrenchment, then both mutations will be unfavorable, since the site is entrenched at its preferred identity in each homolog.

Figure 8 shows the distribution of effects for mutating all sites that differ between BG505 and BF520 to the identity in the other Env. As expected under entrenchment, the average effect of these mutations is deleterious—although there are a substantial number of sites where the mutational flips are not deleterious. We can get some sense of the magnitude of the entrenchment by comparing the effects of the BG505BF520 mutations to the distribution of effects of all possible amino acid mutations (Figure 8). This comparison shows that even unfavorable inter-Env mutational flips are generally more favorable than random amino acid mutations. Therefore, entrenchment occurs for some but not all substitutions that distinguish BG505 and BF520, and the magnitude of entrenchment is less than the effect of a typical random mutation. Entrenchment of substitutions therefore contributes to some of the mutational shifts. But given that many of these shifts occur at sites that do not even differ between the Envs (Figure 7D), entrenchment of substitutions is clearly not the only cause of the shifting amino acid preferences.

Entrenchment of substitutions during Env evolution.

There are 12,521 possible amino acid mutations at the 659 mutagenized sites alignable between BG505 and BF520. The blue densities show the effects of all these mutations to each Env. The orange densities show the effects of just the 92 mutations that convert BG505 to BF520 or vice versa. In the absence of entrenchment, mutating a site in BG505 to its identity in BF520 should have the opposite effect of mutating the site in BF520 to its identity in BG505. In this case, we would expect the BF520BG505 distribution to be the mirror image of the BG505BF520 distribution—and both distributions should be centered around zero if the two Envs are equivalently functional. Instead, mutating a site in either Env to its identity in the other Env tends to be deleterious, indicating that substitutions are often entrenched in the Env in which they have fixed. The effect of a mutation is quantified as the log of the ratio of the site’s preference for the mutant amino acid to the preference for the wild-type amino acid.

https://doi.org/10.7554/eLife.34420.029

Comparing selection in the lab to natural selection

Our experiments measure the effects of mutations on viral growth in a T-cell line in the lab. But HIV actually evolves in humans, where additional selection pressures on Env are undoubtedly present. For instance, antibody pressure might increase the rate of evolution at some sites (Albert et al., 1990; Wei et al., 2003; Richman et al., 2003), whereas pressure to mask certain epitopes (Kwong et al., 2002) might add constraint at other sites. Comparing selection in our experiments to natural selection can identify sites that are under such additional pressures during HIV’s actual evolution in humans.

We determined whether each site in Env evolves faster or slower in nature than expected given three models: that evolution is purely neutral (all nonsynonymous and synonymous mutations have equivalent effects), that sites are under the protein-level constraint measured in our experiments with BG505, or that sites are under the constraint measured with BF520. The first model used a standard dN/dS test (Kosakovsky Pond and Frost, 2005), whereas the other two models are conceptually similar but account for the experimentally measured amino acid preferences as described by (Bloom, 2017). All three models test if individual sites evolve faster or slower than expected, but they ‘expect’ different things: the dN/dS model expects nonsynonymous and synonymous mutations fix at the same rate, while the ExpCM expects the rate at a site to depend on the experimentally measured functional constraints. In all cases, the evidence that a site r evolves differently than expected is statistically summarized by the p-value that ωr is > or < 1. The standard dN/dS model finds hundreds of sites that evolve slower than expected under neutral evolution (Table 1, ωr<1), and only a handful of sites that evolve faster than expected under neutral evolution (Table 1, ωr>1). This finding is unsurprising, since it is well known that Env is under functional constraint. In contrast, ExpCM’s that test the rates of evolution relative to the experimentally measured constraints find far fewer sites that evolve slower than expected, but many more sites that evolve faster (Table 1).

The sites that evolve slower or faster than expected from the experiments are shown in Figure 9A, B, and overlaid on the logoplots in Figure 4 and Figure 5 as the ωr values. The identified sites are similar regardless of whether we use the experiments with BG505 or BF520 (Figure 9C). The reason the results are similar for both experimental datasets is that (as discussed above) the amino acid preferences of most sites are similar in both Envs, suggesting that either dataset provides a reasonable approximation of the site-specific functional constraints across the clade A Envs in Figure 1.

Figure 9 with 3 supplements see all
Sites in Env that evolve faster or slower in nature than expected given the functional constraints measured in the lab.

We calculated the statistical evidence that each site evolves faster (ωr>1) or slower (ωr<1) than expected given the experimentally measured amino acid preferences using the method of Bloom (2017). (A) One monomer of the Env trimer (Stewart-Jones et al., 2016) is colored from blue to white to red based on the strength of evidence that sites evolve slower than expected (blue), as expected (white) or faster than expected (red) given the BG505 experiments. Sites for which we lack ωr estimates are colored black. Sites where the rate of evolution is significantly different than expected at a false discovery rate of 0.05 are shown in spheres. (B) Like (A) but using the data from the BF520 experiments. For both Envs, sites that evolve significantly slower or faster than expected are often on Env’s surface Figure 9—figure supplement 1. (C) The results are similar regardless of whether the BG505 or BF520 experiments are used. Many of the sites of slower-than-expected evolution are asparagines in N-linked glycosylation motifs Figure 9—figure supplement 2. All sites that evolve slower than expected for both experimental datasets are in Figure 9—figure supplement 3. (D) A large cluster of sites that evolve slower than expected is likely involved in Env’s transition between open and closed conformational states. Gray boxes indicate sites that (Ozorowski et al., 2017) proposed form a hydrophobic network that regulates the conformational change; blue boxes and sticks indicate sites that evolve slower than expected. All analyses used the phylogenetic tree in Figure 1. The ωr and Q-values are in Figure 9— source data 1.

https://doi.org/10.7554/eLife.34420.030

What causes some sites to evolve faster or slower in nature than expected from the experiments? The answer in both cases is likely to be immune selection. Most of the sites of faster-than-expected evolution are on the surface of Env (Figure 9A, B and Figure 9—figure supplement 1). Env’s escape from autologous neutralizing antibodies often involves amino acid substitutions in surface-exposed regions (Moore et al., 2009), including at many of the sites that evolve faster than expected. Since our deep mutational scanning did not impose antibody pressure, sites where substitutions are antibody-driven will evolve faster in nature than expected from the experiments.

Interestingly, immune selection also offers a plausible explanation for the sites that evolve slower than expected. In addition to escaping immunity via substitutions at antibody-binding footprints, Env is notorious for employing a range of more general strategies to reduce its susceptibility to antibodies. These strategies include shielding immunogenic regions with glycans (Wei et al., 2003; Stewart-Jones et al., 2016; Gristick et al., 2016) or hiding them by adopting a closed protein conformation (Kwong et al., 2002; Guttman et al., 2015; Ozorowski et al., 2017). Sites that contribute to such general immune-evasion strategies will be under a constraint in nature that is not present in our experiments—and indeed, such sites evolve more slowly than expected from our experiments. For instance, we find very little selection to maintain most glycans in our cell-culture experiments. Of the 21 N-linked glycosylation sites shared between BG505 and BF520, only four are under strong selection to maintain the glycan in our experiments—despite the fact that most are conserved in nature (Figure 9C and Figure 9—figure supplement 2). This finding concords with prior literature suggesting that these glycans are selected primarily for their role in immune evasion (Pugach et al., 2004; Wang et al., 2013; Rathore et al., 2017). Similarly, a network of sites that help regulate Env’s transition between open and closed conformations that have different antibody susceptibilities (Figure 9D) also evolve slower in nature than expected from our experiments. Therefore, we can distinguish evolutionary patterns that are shaped by simple selection for Env function from those that are due to the additional complex pressures imposed during human infections.

Discussion

We have experimentally measured the preference for each amino acid at each site in the ectodomain and transmembrane domain of two Envs under selection for viral growth in cell culture. These amino acid preference maps are generally consistent with prior knowledge about sites that are important for protein properties such as receptor binding or disulfide-mediated stability. However, the main value of these maps comes not from comparing them with prior knowledge, but from the fact that such prior knowledge encompasses just a small fraction of the vast mutational space available to Env. Because Env evolves so rapidly, every study of this protein must be placed in an evolutionary context, and our comprehensive amino acid preference maps potentially enable this in ways that prior piecemeal studies of mutations cannot.

But these maps come with a potentially serious caveat: each one is measured for just a single Env variant. The major question that our study aimed to answer is whether the maps are still useful for evolutionary questions, or whether Env’s amino acid preferences shift so rapidly that each map only applies to the specific HIV strain for which it was measured. This question is reminiscent of one that was grappled with in the early days of protein crystallography, when it first became possible to build maps of a protein’s structure. Because it was not (and is still not) possible to crystallize every variant of a protein, it was necessary to determine whether protein structures could be usefully generalized among homologs. Fortunately for the utility of structural biology, it soon became apparent that closely homologous proteins have similar structures (Chothia and Lesk, 1986; Sander and Schneider, 1991). This rough generalizability of protein structures holds even for a protein as conformationally complex as Env—for although there are many examples of mutations that alter aspects of Env’s conformation and dynamics (Kwong et al., 2000; White et al., 2010; Almond et al., 2010; Davenport et al., 2013), SOSIP trimer structures from diverse Env strains remain highly similar in most respects (Julien et al., 2015; Pugach et al., 2015; Stewart-Jones et al., 2016; Verkerke et al., 2016; Gristick et al., 2016).

Our results show that amino acid preference maps of Env also have a useful level of conservation for many purposes. From a qualitative perspective, the amino acid preferences look mostly similar between BG505 and BF520, and so provide a valuable reference for estimating which mutations are likely to be tolerated at each site in diverse HIV strains. Indeed, we anticipate that the complete maps of mutational effects in Figure 4 and Figure 5 will be useful for future sequence-structure-function studies. From an analytical perspective, a powerful use of our maps is to identify sites that evolve differently in nature than is required by the simple selection for viral growth imposed in our experiments—and the identified sites are largely the same regardless of whether the analysis uses an amino-acid preference map from BG505 or BF520.

Of course, from the perspective of protein evolution, the most interesting sites are the exceptions to the general conservation of amino-acid preferences. Consistent with studies of other proteins (Natarajan et al., 2013; Harms and Thornton, 2014; Doud et al., 2015; Starr et al., 2017), we find a subset of sites that change markedly in which mutations they tolerate. Some shifted sites simply accommodate more amino acids in the more stable BG505 Env—a type of shift that has been well-documented for other proteins (Wang et al., 2002; Bloom et al., 2006; Gong et al., 2013; Kumar et al., 2017). But interestingly, there is no strong trend for shifts to be enhanced at sites that differ between BG505 and BF520. Recent studies of protein evolution have focused on the idea that substitutions become ‘entrenched’ as sites shift to accommodate new amino acids (Pollock et al., 2012; Shah et al., 2015; Bazykin, 2015; Starr et al., 2017). Indeed, a prior protein-wide comparison of amino acid preferences across homologs of influenza nucleoprotein found a significant enrichment of shifts at sites of substitutions (Doud et al., 2015). But although there is some entrenchment of differences between BG505 and BF520, this is not the major factor behind the shifts in amino acid preferences: the majority of sites that have shifted between BG505 and BF520 actually have the same wild-type amino acid in both Envs even though the preferences have shifted. This rather surprising result might be due to Env’s exceptional conformational complexity—mutations can cause long-range alterations in Env’s conformation (Kwong et al., 2000; White et al., 2010; Almond et al., 2010; Davenport et al., 2013), so it seems plausible that they might also shift mutational tolerance at distant sites. Regardless of the exact mechanism, our large-scale datasets of mutational effects in multiple viral strains should be useful for efforts to computationally parameterize ‘fitness landscapes’ of Env (Kouyos et al., 2012; Ferguson et al., 2013; Mann et al., 2014; Barton et al., 2015; Louie et al., 2018).

Our experiments provide highly quantitative data on the mutational tolerance of Env under selection for viral growth in cell culture. These data are amenable to rigorous functional and evolutionary analyses. Here, we have shown how these data can be compared between Envs to identify sites where mutational tolerance shifts with viral genotype, or between experiments and nature to identify sites under different pressure in the lab and in humans. Future experiments that modulate selection pressures in other relevant ways should provide further insight into the forces that drive and constrain HIV’s evolution.

Materials and methods

Creation of codon-mutant libraries

Our codon mutant libraries mutagenized all sites in env to all 64 codons, except that the signal peptide and cytoplasmic tail were not mutagenized. The rationale for excluding these regions is that they are not part of Env’s ectodomain and are prone to mutations that strongly modulate Env’s expression level (Chakrabarti et al., 1989; Yuste et al., 2004; Li et al., 1994).

The codon-mutant libraries were generated using the approach originally described in (Bloom, 2014a), with the modification of (Dingens et al., 2017) to ensure more uniform primer melting temperatures. The computer script used to design the mutagenesis primers (along with some detailed implementation notes) is at https://github.com/jbloomlab/CodonTilingPrimers. For BF520, the three libraries are the same ones described by (Dingens et al., 2017). For BG505, we created three libraries for this study. The wild-type BG505 sequence used for these libraries is in Supplemental file 3. The BG505 mutagenesis primers are in Supplemental file 4.

The end primers for the BG505 mutagenesis were: 5’-tgaaggcaaaactactggtccgtctcgagcagaagacagtggcaatgaga-3’ and 5’-gctacaaatgcatataacagcgtctcattctttccctaacctcaggcca-3’. As with BF520, we cloned the BG505 env libraries into the env locus of the full-length proviral genome of HIV strain Q23 (Poss and Overbaugh, 1999) using the high-efficiency cloning vector described in (Dingens et al., 2017). For this cloning, we digested the cloning vector with BsmBI, and then used PCR to elongate the amplicons to include 30 nucleotides at each end that were identical in sequence to the ends of the BsmBI-digested vector. The primers for this PCR were: 5’-agataggttaattgagagaataagagaaagagcagaagacagtggcaatgagagtgatgg-3’ and 5’-ctcctggtgctgctggaggggcacgtctcattctttccctaacctcaggccatcc-3’. Next, we used NEBuilder HiFi DNA Assembly (NEB, E2621S) to clone the env amplicons into the BsmBI-digested plasmids. We purified the assembled products using Agencourt AMPure XP beads (Beckman Coulter, A63880) using a bead-to-sample ratio of 1.5, and then transformed the purified products into Stellar electrocompetent cells (Takara, 636765). The transformations yielded between 1.5 and 3.6 million unique clones for each of the three replicate libraries, as estimated by plating 1:2000 dilutions of the transformations. We scraped the plated colonies and maxiprepped the plasmid DNA; unlike in (Dingens et al., 2017), we did not include a 4 hr outgrowth step after the scraping step. For the wild-type controls, we maxiprepped three independent cultures of wildtype BG505 env cloned into the same Q23 proviral plasmid. See Figure 2—figure supplement 1 and Figure 3A for information on the average mutation rate in these libraries as estimated by Sanger sequencing and deep sequencing, respectively.

Generation and passaging of viruses

For BG505, we generated mutant virus libraries from the proviral plasmid libraries by transfecting 293 T cells in three 6-well plates (so 18 wells total per library) with a per-well mixture of 2 μg plasmid DNA, 6 μl FuGENE 6 Transfection Reagent (Promega, E269A), and 100 μl DMEM. The 293 T cells were seeded at 5×105 cells/well in D10 media (DMEM supplemented with 10% FBS, 1% 200 mM L-glutamine, and 1% of a solution of 10,000 units/mL penicillin and 10,000 μg/mL streptomycin) the day before transfection, such that they were approximately 50% confluent at the time of transfection. In parallel, we generated wildtype viruses by transfecting one six-well plate of 293 T cells with each wildtype plasmid replicate. At 2 days post-transfection, we harvested the transfection supernatant, passed it through a 0.2 μm filter to remove cells, treated the supernatant with DNAse to digest residual plasmid DNA as in (Haddox et al., 2016), and froze aliquots at −80C. We thawed and titered aliquots using the TZM-bl assay in the presence of 10 μg/mL DEAE-dextran as described in (Dingens et al., 2017).

We conducted the low MOI viral passage illustrated in Figure 2A in SupT1.CCR5 cells (obtained from Dr. James Hoxie; Boyd et al., 2015). The SupT1.CCR5 cells tested negative for mycoplasma. The SupT1.CCR5 cell line was previously created by engineering the parental SupT1 cell line to express CCR5 (Boyd et al., 2015). We used antibody staining followed by flow cytometry to validate that our stock of SupT1.CCR5 cells expressed CCR5, CXCR4, and CD4. There is no validated STR profile for SupT1.CCR5 cells. However, we performed STR profiling on our stock of cells and compared the results to the ATCC SupT1 (ATCC #CRL-1942) reference profile. We found that 11 of 14 alleles plus both amelogenin alleles matched the reference, with no additional mismatched alleles ins the SupT1.CCR5 profile. Given the known instability of lymphoma cell lines (Inoue et al., 2000), this level of identity suggests that the SupT1.CCR5 cells are indeed related to the parental SupT1 cells (Capes-Davis et al., 2013)

During this passage, cells were maintained in R10 media, which has the same composition as the D10 described above, except RPMI-1640 (GE Healthcare Life Sciences, SH30255.01) is used in the place of DMEM. In addition, the media contained 10 μg/mL DEAE-dextran to enhance viral infection. We infected cells with 4 million (for replicate 1) or 5 million (for replicates 2 and 3) TZM-bl infectious units of mutant virus at an MOI of 0.01, with cells at a starting concentration of 1 million cells/mL in vented tissue-culture flasks (Fisher Scientific, 14-826-80). At day one post-infection, we pelleted cells, aspirated the supernatant, and resuspended cell pellets in the same volume of fresh media still including the DEAE-dextran. At 2 days post-infection, we doubled the volume of each culture with fresh media still including DEAE-dextran. At 4 days post-infection, we pelleted cells, passed the viral supernatant through a 0.2 μm filter, concentrated the virus 30 fold using ultracentrifugation as described in (Dingens et al., 2017), and froze aliquots at −80C. In parallel, for each replicate, we also passaged 2×105 (for replicate 1) or 5×105 (for replicates 2 and 3) TZM-bl infectious units of wildtype virus using the same procedure. To obtain final titers for our concentrated virus, we thawed one of the aliquots stored at −80C and titered using the TZM-bl assay in the presence of 10 μg/mL DEAE-dextran.

For the final short-duration infection illustrated in Figure 2A, for each replicate we infected 106 TZM-bl infectious units into 106 SupT1.CCR5 cells in the presence of 100 μg/mL DEAE-dextran (note that this is a 10-fold higher concentration of DEAE-dextran than for the other steps, meaning that the effective MOI of infection is higher if DEAE-dextran has the expected effect of enhancing viral infection). Three hours post-infection, we pelleted the cells and resuspended them in fresh media without any DEAE-dextran. At 12 hr post-infection, we pelleted cells, washed them once with PBS, and then used a miniprep kit to harvest reverse-transcribed unintegrated viral DNA (Haddox et al., 2016).

The generation, passaging and deep sequencing of BF520 was done in a highly similar fashion, except that we only had a single replicate of the wild-type control. Note that the final passaged BF520 mutant libraries analyzed here actually correspond to the ‘no-antibody’ controls described in (Dingens et al., 2017), but that study did not analyze the initial plasmid mutant libraries relative to these passaged viruses, and so was not able to provide measurements of the amino acid preferences.

Illumina deep sequencing

We deep sequenced all of the samples shown in Figure 3A: the plasmid mutant libraries and wildtype plasmid controls, and the cDNA from the final mutant viruses and wildtype virus controls. In order to increase the sequence accuracy, we used a barcoded-subamplicon sequencing strategy. This general strategy was originally applied in the context of deep mutational scanning by Wu et al. (2014), and the specific protocol used in our work is described in Doud and Bloom (2016) (see also https://jbloomlab.github.io/dms-tools2/bcsubamp.html).

The primers used for BG505 are in Supplementary file 5. The primers used for BF520 are in (Dingens et al., 2017). The data generated by the Illumina deep sequencing are on the Sequence Read Archive under the accession numbers provided at the beginning of the Jupyter notebook in Supplementary file 1 and 2.

Analysis of deep-sequencing data

We analyzed the deep-sequencing data using the dms_tools2 software package (Bloom, 2015https://jbloomlab.github.io/dms_tools2/, version 2.2.4). The algorithm that goes from the deep-sequencing counts to the amino acid preferences is that described in (Bloom, 2015) (see also https://jbloomlab.github.io/dms-tools2/prefs.html). A Jupyter notebook that performs the entire analysis including generation of most of the figures in this paper is in Supplementary file 1. An HTML rendering of this notebook is in Supplementary file 2. A repository containing all of this code is also available at https://github.com/jbloomlab/EnvMutationalShiftsPaper (Haddox et al., 2018, copy archived at https://github.com/elifesciences-publications/EnvMutationalShiftsPaper).

The Jupyter notebooks in Supplementary file 1 and 2 also contain numerous plots that summarize relevant aspects of the deep sequencing such as read depth, per-codon mutation frequency, mutation types, etc. Supplementary file 1 also contains text files and CSV files with the numerical values shown in these plots.

Citations are also owed to weblogo (Crooks et al., 2004) and ggseqlogo (Wagih, 2017), which were used in the generation of the logoplots.

Alignments and phylogenetic analyses of Env sequences

A basic description of the process used to generate the clade A sequence alignment in tree-source data 1, the alignment mask in tree-source data 1, and the phylogenetic tree in tree are provided in the legend to that figure. An algorithmic description of how the alignment and tree were generated are in Supplementary file 1 and 2.

For fitting of the phylogenetic substitution models, we used Table 1 (Hilton et al., 2017http://jbloomlab.github.io/phydms/, version 2.2.1) to optimize the substitution model parameters and branch lengths on the fixed tree topology intree. The Goldman-Yang (or YNGKP) model used in Table 1 is the M5 variant described by Yang et al. (2000), with the equilibrium codon frequencies determined empirically using the CF3 × 4 method (Kosakovsky Pond et al., 2010). For the ExpCM shown in Table 1, we extended the models with empirical nucleotide frequencies described in Hilton et al. (2017) to also allow ω to be drawn from discrete gamma-distributed categories exactly as for the M5 model. These ExpCM with gamma-distributed ω were implemented in Table 1 using the equations provided by (Yang, 1994) (see also http://jbloomlab.github.io/phydms/implementation.html#models-with-a-gamma-distributed-model-parameter). The preferences were re-scaled by the stringency parameters in Table 1 as described in Hilton et al. (2017). For both the M5 model and the ExpCM with a gamma-distributed ω, we used four categories for the discretized gamma distribution.

Table 1—source data 1 shows the results for a wider set of models than those used in Table 1. These include the M0 model of (Yang et al., 2000), ExpCM without a gamma-distributed ω, and ExpCM in which the amino acid preferences are averaged across sites as a control to ensure that the improved performance of these models is due to their site-specificity. Note how for these Env alignments, using a gamma-distributed ω is very important in order for the ExpCMs to outperform the M5 model—we suspect this is because there are many sites of strong diversifying selection.

For detection of sites with faster or slower than expected evolution, we used the approach in (Bloom, 2017), which is exactly modeled on the FEL approach of (Kosakovsky Pond and Frost, 2005) but extended to ExpCM. This approach estimates a p-value that ωr is not equal to one for each site r using a likelihood-ratio test. The actual point estimates of ωr are unreliable for individual sites due to the limited number of observations, so we report the p-value that ωr is not equal to one, which is a better indication of the strength of the statistical evidence for faster or slower than expected evolution (Kosakovsky Pond and Frost, 2005; Murrell et al., 2012). For the Q-values and false discovery rate testing, we considered the tests for ωr>1 and ωr<1 separately.

Supplementary file 1 and 2 contains the code that runs Table 1 to reproduce all of these analyses.

Re-scaling the preferences

The amino acid preferences that are directly extracted from the deep sequencing data essentially give the enrichment/depletion of each mutation, normalized to sum to one at each site (Doud et al., 2015, https://jbloomlab.github.io/dms_tools2/prefs.html). However, the extent that any mutation is enriched or depleted is a combination of two factors: the inherent effect of that mutation, and the ‘stringency’ of the experimental selection. For instance, if the selection is weak, then deleterious mutations will only be slightly depleted; conversely, if selection is strong, then deleterious mutations will be greatly depleted. The fact that the preferences depend on the stringency of the experimental selection is important if we want to compare results between Envs. The reason is that our goal is to identify differences in the inherent effects of mutations between Envs, not simply find differences due to variation in experimental stringency. Of course, we have done our best to perform the experiments for BG505 and BF520 equivalently, but because these are different viruses with different growth rates, it is impossible to exactly match the experimental stringencies. This can be seen in Figure 3—source data 1, which shows that stop codons were more depleted for BG505 than BF520, indicating that selection in our experiments was more stringent for BG505.

How should we best re-scale the preferences? Raising them to a power is a sensible approach. To see why, imagine a mutation that is depleted 3-fold after 2 rounds of viral growth. If our experiment instead allowed 22=4 rounds of viral growth, then the mutation would be depleted 32=9 -fold. More generally, if a mutation is enriched in frequency by ϕ-fold after n rounds of viral growth, then it will be enriched in frequency by ϕβ-fold after β×n rounds of viral growth. Since the amino acid preferences are conceptually equivalent to the re-normalized enrichments of mutations (Bloom, 2015) https://jbloomlab.github.io/dms-tools2/prefs.html, it therefore makes sense that the re-scaled preference πr,a for amino acid a at site r should be related to the directly measured preference π^r,a by πr,a(π^r,a)β. And indeed, this is exactly the re-scaling scheme described in (Hilton et al., 2017) that we use to re-scale our preferences for BG505 and BF520.

The last point is how to choose the re-scaling parameter β for each Env. It turns out that the features that we have described above for our experiments are also a feature of natural evolution: the expected frequency of a substitution during evolution depends not only on the inherent fitness effect of that mutation, but also on the effective population size, which is conceptually somewhat similar to the stringency of selection. It turns out that in a mutation-selection phylogenetic model of evolution, if the amino acid preferences are taken to represent the ‘fitness effects’ of mutations, then the exponential scaling parameter β is proportional to the effective population size (Halpern and Bruno, 1998; McCandlish and Stoltzfus, 2014; Bloom, 2014b). Therefore, fitting the β parameter using a phylogenetic approach enables standardization of the preferences for the two Envs, and re-scales the preferences so that they best match with the actual stringency of selection observed in nature (Hilton et al., 2017).

Note that in practice this re-scaling scheme is roughly equivalent to a more heuristic approach that has been used by (Gray et al., 2017) and others. In this heuristic approach, the log-transformed enrichment ratios from different experiments are adjusted so that the distributions have equal spreads. Since multiplying log-transformed enrichment ratios is equivalent to exponentiating amino acid preferences, these two re-scaling procedures apply the same mathematical transformation.

Identifying sites of shifted amino acid preference

When identifying shifts in amino acid preferences between the two Envs, we needed a way to quantify differences between the Envs while accounting for the fact that our measurements are noisy. The approach we use is based closely on that of Doud et al. (2015) and is illustrated graphically in Figure 6A. The RMSDcorrected value is our measure of the magnitude of the shift. Figure 6A, its legend, and the associated text completely explains these calculations with the following exception: they do not detail how the ‘distance’ between any two preference measurements was calculated. The distance between preferences at each site was simply defined as half of the sum of absolute value of the difference between preferences for each amino acid. Specifically, for a given site r, let πr,ai be the re-scaled preference for amino acid a in homolog i (e.g. BG505) and let πr,aj be the re-scaled preference for that same amino acid in homolog j (e.g. BF520). Then the distance between the homologs at this site is simply Dri,j=12a|πr,aiπr,aj|. The factor of 12 is used so that the maximum distance will always fall between zero and one.

Analysis of entrenchment

For the analysis in Figuer 8, the results are presented in terms of the mutational effects rather than the amino acid preferences. If πr,a is the preference of site r for amino acid a and πr,a is the preference for amino acid a (both re-scaled by the stringency parameters in Table 1), then the estimated effect of the mutation from a to a is simply log(πr,aπr,a).

Data and code availability

All code and input data required to reproduce all analyses in this paper are in Supplementary file 1 (see also Supplementary file 2). A repository containing all of this code is also available at https://github.com/jbloomlab/EnvMutationalShiftsPaper (Haddox et al., 2018; copy archived at https://github.com/elifesciences-publications/EnvMutationalShiftsPaper). The deep sequencing data are on the Sequence Read Archive with the accession numbers listed in Supplementary file 1 and 2.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
    Controlling the false discovery rate: a practical and powerful approach to multiple testing
    1. Y Benjamini
    2. Y Hochberg
    (1995)
    Journal of the Royal Statistical Society, Series B pp. 289–300.
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
    The cytoplasmic domain of simian immunodeficiency virus transmembrane protein modulates infectivity
    1. L Chakrabarti
    2. M Emerman
    3. P Tiollais
    4. P Sonigo
    (1989)
    Journal of Virology 63:4395–4403.
  18. 18
  19. 19
  20. 20
    The relation between the divergence of sequence and structure in proteins
    1. C Chothia
    2. AM Lesk
    (1986)
    The EMBO Journal 5:823.
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
    Numbering positions in HIV relative to HXB2CG
    1. B Korber
    2. BT Foley
    3. C Kuiken
    4. SK Pillai
    5. JG Sodroski
    (1998)
    Human Retroviruses and AIDS 3:102–111.
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
    Initial sequencing and comparative analysis of the mouse genome
    1. RH Waterston
    2. K Lindblad-Toh
    3. E Birney
    4. J Rogers
    5. JF Abril
    6. P Agarwal
    7. R Agarwala
    8. R Ainscough
    9. M Alexandersson
    10. P An
    11. SE Antonarakis
    12. J Attwood
    13. R Baertsch
    14. J Bailey
    15. K Barlow
    16. S Beck
    17. E Berry
    18. B Birren
    19. T Bloom
    20. P Bork
    21. M Botcherby
    22. N Bray
    23. MR Brent
    24. DG Brown
    25. SD Brown
    26. C Bult
    27. J Burton
    28. J Butler
    29. RD Campbell
    30. P Carninci
    31. S Cawley
    32. F Chiaromonte
    33. AT Chinwalla
    34. DM Church
    35. M Clamp
    36. C Clee
    37. FS Collins
    38. LL Cook
    39. RR Copley
    40. A Coulson
    41. O Couronne
    42. J Cuff
    43. V Curwen
    44. T Cutts
    45. M Daly
    46. R David
    47. J Davies
    48. KD Delehaunty
    49. J Deri
    50. ET Dermitzakis
    51. C Dewey
    52. NJ Dickens
    53. M Diekhans
    54. S Dodge
    55. I Dubchak
    56. DM Dunn
    57. SR Eddy
    58. L Elnitski
    59. RD Emes
    60. P Eswara
    61. E Eyras
    62. A Felsenfeld
    63. GA Fewell
    64. P Flicek
    65. K Foley
    66. WN Frankel
    67. LA Fulton
    68. RS Fulton
    69. TS Furey
    70. D Gage
    71. RA Gibbs
    72. G Glusman
    73. S Gnerre
    74. N Goldman
    75. L Goodstadt
    76. D Grafham
    77. TA Graves
    78. ED Green
    79. S Gregory
    80. R Guigó
    81. M Guyer
    82. RC Hardison
    83. D Haussler
    84. Y Hayashizaki
    85. LW Hillier
    86. A Hinrichs
    87. W Hlavina
    88. T Holzer
    89. F Hsu
    90. A Hua
    91. T Hubbard
    92. A Hunt
    93. I Jackson
    94. DB Jaffe
    95. LS Johnson
    96. M Jones
    97. TA Jones
    98. A Joy
    99. M Kamal
    100. EK Karlsson
    101. D Karolchik
    102. A Kasprzyk
    103. J Kawai
    104. E Keibler
    105. C Kells
    106. WJ Kent
    107. A Kirby
    108. DL Kolbe
    109. I Korf
    110. RS Kucherlapati
    111. EJ Kulbokas
    112. D Kulp
    113. T Landers
    114. JP Leger
    115. S Leonard
    116. I Letunic
    117. R Levine
    118. J Li
    119. M Li
    120. C Lloyd
    121. S Lucas
    122. B Ma
    123. DR Maglott
    124. ER Mardis
    125. L Matthews
    126. E Mauceli
    127. JH Mayer
    128. M McCarthy
    129. WR McCombie
    130. S McLaren
    131. K McLay
    132. JD McPherson
    133. J Meldrim
    134. B Meredith
    135. JP Mesirov
    136. W Miller
    137. TL Miner
    138. E Mongin
    139. KT Montgomery
    140. M Morgan
    141. R Mott
    142. JC Mullikin
    143. DM Muzny
    144. WE Nash
    145. JO Nelson
    146. MN Nhan
    147. R Nicol
    148. Z Ning
    149. C Nusbaum
    150. MJ O'Connor
    151. Y Okazaki
    152. K Oliver
    153. E Overton-Larty
    154. L Pachter
    155. G Parra
    156. KH Pepin
    157. J Peterson
    158. P Pevzner
    159. R Plumb
    160. CS Pohl
    161. A Poliakov
    162. TC Ponce
    163. CP Ponting
    164. S Potter
    165. M Quail
    166. A Reymond
    167. BA Roe
    168. KM Roskin
    169. EM Rubin
    170. AG Rust
    171. R Santos
    172. V Sapojnikov
    173. B Schultz
    174. J Schultz
    175. MS Schwartz
    176. S Schwartz
    177. C Scott
    178. S Seaman
    179. S Searle
    180. T Sharpe
    181. A Sheridan
    182. R Shownkeen
    183. S Sims
    184. JB Singer
    185. G Slater
    186. A Smit
    187. DR Smith
    188. B Spencer
    189. A Stabenau
    190. N Stange-Thomann
    191. C Sugnet
    192. M Suyama
    193. G Tesler
    194. J Thompson
    195. D Torrents
    196. E Trevaskis
    197. J Tromp
    198. C Ucla
    199. A Ureta-Vidal
    200. JP Vinson
    201. AC Von Niederhausern
    202. CM Wade
    203. M Wall
    204. RJ Weber
    205. RB Weiss
    206. MC Wendl
    207. AP West
    208. K Wetterstrand
    209. R Wheeler
    210. S Whelan
    211. J Wierzbowski
    212. D Willey
    213. S Williams
    214. RK Wilson
    215. E Winter
    216. KC Worley
    217. D Wyman
    218. S Yang
    219. SP Yang
    220. EM Zdobnov
    221. MC Zody
    222. ES Lander
    223. Mouse Genome Sequencing Consortium
    (2002)
    Nature 420:520–562.
    https://doi.org/10.1038/nature01262
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
    Identification of individual human immunodeficiency virus type 1 gp120 amino acids important for CD4 receptor binding
    1. U Olshevsky
    2. E Helseth
    3. C Furman
    4. J Li
    5. W Haseltine
    6. J Sodroski
    (1990)
    Journal of Virology 64:5701–5707.
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
    Variants from the diverse virus population identified at seroconversion of a clade A human immunodeficiency virus type 1-infected woman have distinct biological properties
    1. M Poss
    2. J Overbaugh
    (1999)
    Journal of Virology 73:5255–5264.
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83
  84. 84
  85. 85
  86. 86
  87. 87
  88. 88
  89. 89
  90. 90
  91. 91
  92. 92
  93. 93
  94. 94
  95. 95
  96. 96
  97. 97
  98. 98
  99. 99
  100. 100
  101. 101
  102. 102
  103. 103
  104. 104
  105. 105
  106. 106
  107. 107
  108. 108
  109. 109
  110. 110
  111. 111
  112. 112
    Codon-substitution models for heterogeneous selection pressure at amino acid sites
    1. Z Yang
    2. R Nielsen
    3. N Goldman
    4. AM Pedersen
    (2000)
    Genetics 155:431–449.
  113. 113
  114. 114
  115. 115
  116. 116

Decision letter

  1. Arup K Chakraborty
    Reviewing Editor; Massachusetts Institute of Technology, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Mapping mutational effects along the evolutionary landscape of HIV envelope" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by Arup Chakraborty as the Reviewing Editor and Senior Editor. The following individuals involved in review of your submission has agreed to reveal their identity: John P Barton (Reviewer #1); Peter Kim (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Information on the possible amino acid trajectories of the rapidly evolving HIV envelope protein has long been limited. In this manuscript, you extend previous work from your lab analyzing the mutational landscape of the HIV-1 Env gene during in vitro viral replication. The work focuses on two clade A transmitter founder viruses, BG505 and BF520, which differ by more than 100 mutations. A deep mutational scanning approach is used to identify the permitted mutations at each residue of Env. Of particular interest is how amino acid preferences differ for these two proteins. A set of 30 sites with clear shifts in amino acid preference are identified. These data are compared to known mutational constraints on Env and to one another, and sources for the shifts in mutant affects are proposed by eliminating many scenarios. The results offer a body of knowledge useful for identifying i) residues and networks that are critical for basic function and ii) hot spots for evolvability or immune pressure. The analysis has been performed to a high standard (though some clarifications are necessary), and the openness of the data is a great service to the community. The paper represents a significant contribution towards understanding the structure-function relationship of HIV envelope. Therefore, the paper is likely to be accepted for publication if the following points are adequately addressed.

Essential revisions:

1) How are mutations that strongly modulate Env's expression level taken into account?

2) Only the native, pre-fusion Env trimer is used to calculate proximity of substituted and shifted residues. It is possible that substitution-shift colocalization may be detected in other conformations of Env as it undergoes structural rearrangement to mediate fusion, as opposed to acting through longer-range allostery in the pre-fusion conformation. Another well-characterized structure that should be considered is the post-fusion six-helix bundle of gp41.

3) Overall, it is found that the shifted sites tend to cluster together, but they are not necessarily the same as the sites with substitutions between the two Envs or ones that contact them. This analysis is not pursued in great detail, however, leaving the reasons why shifted sites tend to cluster together unexplained. Given that the authors also observe some evidence of entrenchment of substitutions between the Envs, is it possible that interactions between the shifted sites themselves are important? This question could be explored through a more detailed examination of the structure.

https://doi.org/10.7554/eLife.34420.046

Author response

Essential revisions:

1) How are mutations that strongly modulate Env's expression level taken into account?

This is an important point. Our experiments directly measure how mutations affect viral growth. This approach has many advantages over simple surface display approaches, since it is more relevant to the true function of Env. However, a caveat is that viral growth is a convolution of Env function and expression, so our experiments cannot separate these two phenotypes.

We have taken steps to limit simply finding mutations affect expression. We have done this by not mutagenizing the signal peptide and cytoplasmic tail, which are the most common sites of mutations that affect expression levels.

We have added text to make this point clear, and emphasize the associated caveats. When we first describe the deep mutational scanning approach, we now say:

“However, in our experiments we mutagenized only the ectodomain and transmembrane domain of Env, and excluded the signal peptide and cytoplasmic tail. The reason is that we measure how Env mutations affect viral growth, which is influenced both by the functionality of Env protein molecules and their expression level. Mutations in the signal peptide and cytoplasmic tail commonly affect Env expression level (Chakrabarti et al., 1989; Yuste et al., 2004; Li et al., 1994), so we excluded these regions with the goal of reducing the degree to which we simply identified mutations that affected Env expression.”

Then when we describe the results, we clearly emphasize how expression level (like Env protein stability, which was already discussed) could have general effects on the results:

“Differences in Env’s expression level might also contribute to a general broadening or narrowing of tolerance to subsequent mutations. Because our experiments select for viral growth – which is affected by both Env function and expression – it is possible that some of the shifts are due to epistatic mutational effects on expression rather than function.”

2) Only the native, pre-fusion Env trimer is used to calculate proximity of substituted and shifted residues. It is possible that substitution-shift colocalization may be detected in other conformations of Env as it undergoes structural rearrangement to mediate fusion, as opposed to acting through longer-range allostery in the pre-fusion conformation. Another well-characterized structure that should be considered is the post-fusion six-helix bundle of gp41.

This is an excellent suggestion. In the initial paper, we only analyzed the closed prefusion conformation of the Env trimer. We now analyze substitution-shift co-localization in two other conformations: the open CD4-bound conformation and the post-fusion six-helix bundle. In the end, our basic conclusions remain the same, adding greater support to the hypothesis that epistasis is acting via long-range interactions.

First, we investigated the open CD4-bound conformation. Figure 7D statistically tests the hypothesis that shifts are larger at sites that have substituted or sites that contacts substitutions. Initially, this figure only showed the results of this test for the closed prefusion conformation. Now, it shows the results for both the closed and open conformations. Both conformations gave the same general result: there is no significant trend for shifted sites to contact sites of substitutions. We also considered both structures at once, classifying sites as contacting substitutions if they do so in at least one structure. Even then, we did not find an association (see the new right-most panel in Figure 7D). Thus, our investigation of the open conformation re-enforced our initial conclusion that many of the shifts are likely due to long-range, rather than short-range, epistatic interactions.

We also repeated the analyses in Figure 7B and Figure 7C with the open conformation. Our initial conclusions about the solvent accessibility (Figure 7B) and clustering of shifted sites (Figure 7B) also held for the open conformation. Next, we investigated substitution-shift co-localization in the post-fusion six-helix bundle. Only ∼80 residues per Env monomer are resolved in available structures of the six-helix bundle. The structure with largest number of shifted sites (PDB: 1ENV) only contains four shifted sites and five substituted sites. These numbers are too low for the statistical tests used in Figure 7. But we have added a new figure (Figure 7—figure supplement 2) that qualitatively analyzes substitution-shift co-localization in the post-fusion six-helix bundle. Three of the four shifted sites cluster at one end of the helix, along with one of the five substitutions. However, this cluster is also present in the closed and open states as well (Figure 7—figure supplement 2). Thus, the post-fusion six-helix bundle did not reveal any new patterns of co-localization that were not also present in the other states.

We now describe the results of the analyses in the first two paragraphs of the Results subsection “Structural and evolutionary properties of shifted sites”, interleaved with the original text in these paragraphs.

In the process of making these revisions, we discovered a minor bug that affected our calculations of residue-residue distances. We fixed that bug in the revised version. None of the main results change, but some of the quantitative values in the box plots and the corresponding P-values have changed slightly.

3) Overall, it is found that the shifted sites tend to cluster together, but they are not necessarily the same as the sites with substitutions between the two Envs or ones that contact them. This analysis is not pursued in great detail, however, leaving the reasons why shifted sites tend to cluster together unexplained. Given that the authors also observe some evidence of entrenchment of substitutions between the Envs, is it possible that interactions between the shifted sites themselves are important? This question could be explored through a more detailed examination of the structure.

To address this comment, we examined the clusters of shifted sites in greater detail. Two clusters occur within conformationally dynamic regions of Env. One hypothesis as to why shifted sites cluster in these regions is that their conformationally dynamic nature allows long-range epistatic interactions to be propagated between shifted and substituted sites that are distant from one another.

We have added a new figure supplement (Figure 7—figure supplement 3) that provides a fine-grained structural view of each of the clusters. We have also added an entire paragraph detailing these observations at the end of the section in the Results called “Structural and evolutionary properties of shifted sites”. Overall, these additions strengthen the manuscript by providing a better structural intuition for how long-range epistatic interactions might occur in Env.

https://doi.org/10.7554/eLife.34420.047

Article and author information

Author details

  1. Hugh K Haddox

    1. Basic Sciences Division and Computational Biology Program, Fred Hutchinson Cancer Research Center, Seattle, United States
    2. Molecular and Cellular Biology PhD program, University of Washington, Seattle, United States
    Contribution
    Conceptualization, Formal analysis, Investigation, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Adam S Dingens
    Competing interests
    No competing interests declared
  2. Adam S Dingens

    1. Basic Sciences Division and Computational Biology Program, Fred Hutchinson Cancer Research Center, Seattle, United States
    2. Molecular and Cellular Biology PhD program, University of Washington, Seattle, United States
    Contribution
    Conceptualization, Investigation, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Hugh K Haddox
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9603-9409
  3. Sarah K Hilton

    1. Basic Sciences Division and Computational Biology Program, Fred Hutchinson Cancer Research Center, Seattle, United States
    2. Department of Genome Sciences, University of Washington, Seattle, United States
    Contribution
    Software, Formal analysis, Writing—review and editing
    Competing interests
    No competing interests declared
  4. Julie Overbaugh

    1. Human Biology Division, Fred Hutchinson Cancer Research Center, Seattle, United States
    2. Epidemiology Program, Fred Hutchinson Cancer Research Center, Seattle, United States
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Writing—review and editing
    Competing interests
    No competing interests declared
  5. Jesse D Bloom

    1. Basic Sciences Division and Computational Biology Program, Fred Hutchinson Cancer Research Center, Seattle, United States
    2. Department of Genome Sciences, University of Washington, Seattle, United States
    Contribution
    Conceptualization, Software, Supervision, Funding acquisition, Investigation, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    jbloom@fredhutch.org
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1267-3408

Funding

National Institutes of Health (R01-AI127893)

  • Jesse D Bloom

National Science Foundation (DGE-1256082)

  • Adam S Dingens

Howard Hughes Medical Institute (Faculty Scholar Grant)

  • Jesse D Bloom

Simons Foundation (Faculty Scholar Grant)

  • Jesse D Bloom

Collaboration for AIDS Vaccine Discovery (OPP1111923)

  • Jesse D Bloom

National Institutes of Health (DP1-DA039543)

  • Julie Overbaugh

National Institutes of Health (T32GM007270)

  • Hugh K Haddox

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Michael Doud and Orr Ashenberg for computer code that formed the basis for some of the analyses. We thank Andrew Ward for pointing out to us that some of the sites with slower-than-expected rates of evolution are involved in Env’s conformational changes upon receptor binding. We thank Kelly Lee for helpful input about the relative stabilities of the BG505 and BF520 Envs. We thank the Fred Hutch Genomics Core for performing the Illumina deep sequencing.

Reviewing Editor

  1. Arup K Chakraborty, Massachusetts Institute of Technology, United States

Publication history

  1. Received: December 17, 2017
  2. Accepted: March 15, 2018
  3. Accepted Manuscript published: March 28, 2018 (version 1)
  4. Version of Record published: April 20, 2018 (version 2)
  5. Version of Record updated: May 8, 2018 (version 3)

Copyright

© 2018, Haddox et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,600
    Page views
  • 257
    Downloads
  • 4
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Evolutionary Biology
    2. Genetics and Genomics
    Cristian Cañestro, Vittoria Roncalli
    Insight
    1. Evolutionary Biology
    2. Genetics and Genomics
    David Jebb, Michael Hiller
    Short Report