Emerging Patterns in HIV-1 gp120 Variable Domains in Anatomical Tissues in the Absence of a Plasma Viral Load

Hiv undergoes extensive diversification within an individual over the course of the infection. The HIV envelope (Env) protein, which contains two subunits gp120 and gp41, exists as a trimer on the surface of the virus. This trimeric structure serves as the point of contact between the virus and the cellular targets of infection: T-cells, macrophages (and their precursors), and dendritic cells.

As a major factor in cell entry, Env is the focus of multiple strategies for drug/vaccine development. The gp120 subunit contains five variable domains (V1–V5), of which three are remarkably elastic (V1, V2, and V4), displaying many nucleotide insertions and deletions (indels) in both intra- and intersubject sequence populations. These indels often result in the gain or loss of potential N-linked glycosylation sites (N-sites). The network of N-sites on gp120 serves as targets for broadly neutralizing antibodies¹ and contributes to immune evasion and escape during HIV-1 infection.²

Understanding the specific selective pressures on the virus in tissue compartments is of critical importance, as increasing evidence suggests that tissues may be a major source of viral rebound after cessation of combined antiretroviral therapy (cART). Although analysis of gp120 variable domains is hampered by extensive indels that disrupt a homologous alignment, a number of interesting observations suggest that patterns of variability reflect changing selective pressures. For example, gp120 tends to be shorter and contain fewer N-sites in early infection³ and may be preferentially transmitted.⁴ Furthermore, cloned early viruses of shorter length replicated efficiently in activated human CD4⁺T cells but not in cultured monocyte-derived macrophages,⁵ suggesting unique replication/selection pressures in tissue-based immune cells.

The flexibility of gp120 variable domains has rarely been assessed in anatomical tissues. In this study, we analyzed gp120 variable domain sequences for length, glycosylation, charge, and structure from anatomical tissues of 7 subjects who were on cART until death with no measurable plasma viral load (denoted NPVL subjects) and 15 subjects who had a detectable plasma viral load at death (denoted PVL subjects). All subjects in the study were infected with HIV Subtype B.

To generate context for length variation found in our study, we downloaded and assessed the length distributions of V1 (n = 5565), V2 (n = 5964), and V4 (n = 8828) Subtype B sequences deposited in the Los Alamos HIV Database (www.hiv.lanl.gov) using the “one sequence per subject” option. While the bulk of the Los Alamos variable regions were restricted to a limited length range (V1 = 24–34 amino acids, V2 = 39–43 amino acids and V4 = 30–33 amino acids), the overall variance for V1, V2, and V4 length was considerable with lengths of 4–58, 23–75, and 6–54 amino acids for each domain, respectively, thus demonstrating the remarkable elasticity of gp120 in the overall population (Fig. 1).

NPVL and PVL subject and sequence information is provided in Table 1. All sequence data were previously published and available in GenBank (JN786871-JN786685, HM001362-HM002482, KU708874-KU70932, KY270561-KY270659).

NPVL participants were enrolled at the National Neurological AIDS Bank (NNAB), a member of the National NeuroAIDS Tissue Consortium (NNTC)⁶ and were on cART until death⁷; it is difficult to obtain such tissues, as the majority of subjects with terminal disease stop taking cART and thus demonstrate a degree of viral rebound. Interestingly, out of the 120 diverse tissues examined in this NPVL cohort, only 24 tissues contained amplifiable HIV env-nef, demonstrating the nonubiquitous nature of HIV in tissues during cART, as well as the ability of cART to reduce the anatomical reservoir, but not eliminate the virus completely.⁸ By contrast, all 76 tissues from the PVL group contained readily amplifiable HIV. The PVL study group comprised autopsy tissue from patients who died with detectable plasma viral load and AIDS-associated disease pathologies.^9–12

All DNA sequences in the study were translated into their coded amino acids; and only those with open reading frames were retained for further analysis. For all sequences, variable domain characteristics (length, number of N-sites, and charge) were calculated using the “Variable Domain Characteristics” tool at the HIV Sequence Database at Los Alamos with additional calculations performed in Excel (ver. 15.20).

Box and whisker plots appeared to show a trend in different length distributions between PVL and NPVL categories (Fig. 2). Therefore, we used a two-tailed t-test (MacWizard ver 1.9.18) to determine whether characteristics of the V1–V5 regions were statistically significantly different between categories at a 95% confidence level. As the number of sequences generated for each subject were uneven, and to avoid biasing the results, in each analysis, the data sets were reduced so that only one sequence of each representative length for each subject was retained for statistical analysis.

There was a significant difference (p < .01) in the distributions of length in V1 between NPVL (mean = 30.93 aa's) and PVL (mean = 28.49 aa's) sequences. Furthermore, the range and distribution of amino acid lengths appeared relatively similar among most of the 24 NPVL tissues, with 43% of tissue's sequence populations containing lengths above and below the min/max quartiles of those in the Los Alamos database shown in Figure 1. With exceptions (notably in subjects HCK1 and HCK2), PVL tissues from specific subjects tended to demonstrate less diversity in V1 length, with the V1 sequences from subject's different tissues usually representing either longer or shorter V1s than those in Los Alamos, and only 8% of the PVL tissue's sequence populations spanned both the min/max quartiles observed in the NPVL subject's V1 sequence populations.

In V2, sequences from NPVL subjects tended to be longer in comparison to PVL subjects. Seventy eight percent of the NPVL tissues from 6/7 subjects harbored sequences above the upper quartile identified in the Los Alamos V2 data, whereas only 23% of tissue sequence populations, representing 4/15 PVL subjects, were above the upper quartile in the Los Alamos V2 data. There was no statistical significance to the patterns observed in V3, V4, or V5 domain sequences (p > .01)(data not shown).

The length of Env variable domains is known to be highly correlated with an increase in N-sites. In the current study, no statistical difference was observed when comparing the mean number of N-sites between NPVL and PVL datasets in the V1 region (p > .01; 2.492 aa's and 2.443 aa's, respectively). In contrast, the mean number of N-sites for the NPVL V2 region was significantly greater than the PVL V2 region (p < .01; 3.676 aa's and 2.027 aa's, respectively).

We then assessed variable domains for charge differences by comparing the overall charge of each variable domain and the frequency of positively charged amino acid residues (K and R) and negatively charged amino acid residues (D and E). In V1, NPVL sequences showed a significant increase in the number of charged amino acid residues in comparison to PVL sequences (p < .01), whereas no significant difference in the number of charge residues was observed between the groups in the V2 region.

In summary, while previous sequence analysis studies of combined V1–V2 loops noted that longer sequences were associated with an increase in N-sites, in the current analysis, we uncoupled these domains and determined that longer V1 sequences were correlated to an increase in the number of charged residues (both positively and negatively charged), whereas increased glycosylation was associated with longer V2 domains.

We were concerned that the sequencing approach could have impacted our findings. Early studies of HIV primarily used a bulk polymerase chain reaction (bPCR) amplification, wherein many targets were amplified together in the same tube before sequencing. While this method could result in resampling and distorted proportions of viral variants, most studies using bPCR reported highly variable sequence populations. Later, especially when HIV researchers started studying subjects with a very low or undetectable viral load, single-genome sequencing (SGS), which uses endpoint dilution to amplify from only one template in a well, became a more common approach to sequence HIV.

To test to see if these two approaches resulted in statistically different length populations in any variable domains studied, we pooled all bPCR-produced data and all SGS-produced data from the PVL sequence population and performed a two-tailed t-test. We found no statistical differences between bPCR and SGS lengths in any variable domain sequence population (p > .01).

We next generated an alignment to observe the amino acid variation calculated for V1 and V2. Only unique sequences were used in each alignment. Because automated alignment programs are unable to handle the diversity in V1 and V2 domains, a manual alignment procedure was developed based on anchoring relatively conserved regions of amino acids in each domain as follows: V1 begins with a highly conserved “CTD” motif, followed by a glycosylated region. While this glycosylated region may contain some charged residues, it is heavily populated with small polar residues, such as those that form putative glycosylation sites, “N,” “T,” and “S.” Toward the 3′ end of V1, there exists a series of charged residues interspersed with miscellaneous and often hydrophobic amino acids. In the alignment procedure, a gap is placed between these two regions following these guidelines: all glycosylation sites must be shifted left of the gap. The first appearance of a charged residue or a methionine that does not inhabit the middle position of an N-site indicates the site where the rest of the amino acids are shifted right.

Very few exceptions were made to these rules; however, <10 sequences contained charged residues and no glycosylation sites, so the insertion site was chosen before the “GEIK” motif at the 3′ end. Alternatively, <10 sequences contained no charged residues between the glycosylation region and the GEIK motif, so the same rule was used.

V2 contains two somewhat conserved domains, also at the 5′ and 3′ ends that were used to anchor the alignment. At the 5′ end, there is rarely any length variation in the first 28 positions, despite considerable amino acid variation. At the 3′ end, a relatively conserved “YRLISC” motif is observed. The length variable region is similar to V1 in that it comprises many small polar residues that generate putative N-sites. Therefore, for V2, a gap was inserted before the YRLISC motif, and all residues for each sequence were shifted left. Spaces were inserted in the V1 and V2 PVL and NPVL data alignments so that the two alignments (PVL and NPVL) correspond to each other. Alignments and graphics shown in Figure 3 were performed using Geneious software (version 9.1.3).

The SeqLogo consensus graphs (Fig. 3), where the height of the amino acid at each site is proportional to the total numbers of identified amino acids in that column, demonstrate the variability observed within each sequence population. The V1 alignment mirrors results observed in the glycosylation and charge analysis: while both groups have N-sites occurring primarily in the first 20 amino acid residues in the flexible region of glycosylation, the trend of longer sequences in the NPVL group is often due to an abundance of charged positions at the 3′ end of the V1 domain. Also, notable is that the amino acid substitutions in the 3′ charged region can result in considerable differences in hydrophobicity and isoelectric point in the V1 loop.

The V2 SeqLogo graph also clearly demonstrates that the primary difference between NPVL and PVL sequence populations lies in length and in residues that contribute to N-sites (N, T, and S). Despite considerable amino acid variation in both V2 sequence populations, no obvious amino acid substitutions occurred that impacted either the hydrophobicity or isoelectric point along the V2 loop.

Functionally, V3 remains sequestered beneath V1–V2 until a target cell is reached, triggering a conformational rearrangement of V1–V2 to an open state, allowing for exposure of V3.¹³ A deeper, highly charged, and glycosylated V3 pocket may mask epitopes recognized by neutralizing antibodies; however, these same variable loop insertions may hinder the conformational process required for CD4 binding.¹⁴ Thus, the modifications observed between long and short V1–V2 domains may have a considerable impact on the function of the envelope.

To observe how length variation in V1 and V2 impacts the structure of the gp120 binding pocket, two sequences representing long V1–V2 domains (GenBank accession numbers KU708975, KU708891) and two sequences representing short V1–V2 sequences (GenBank accession numbers HM002301, HM001816) were chosen for representation as 3D protein structures (Figs. 4 and 5). Structures were modeled using the I-TASSER web server^15,16 and aligned and viewed in MacPyMOL (v1.7.2.1). The top protein threading models were structure 3j70D or 5FYjG from the Protein Data Base (www.rcsb.org/pdb/).

Each pose had a C-score and TM-score within the range of acceptable topologies. In Figure 4, we focus only on the structure of V1–V2. Note that while the overall charge difference between the two sequences is small (−3), the long sequence has twice as many charged residues as the shorter sequence. The structural comparisons clearly demonstrate two other dramatic differences between long and short structures: (1) long sequences have a much deeper binding pocket than short sequences (Fig. 4E), and (2) long sequences have an abundance of charges (contributed by V1) and N-sites (contributed by V2) within the binding pocket (Fig. 4F).

In Figure 5, we characterize two other examples of long and short sequences, this time using the entire V1–V5 structure. Here, the length difference is less than in the structures shown in Figure 4; however, there is still a dramatic contrast between long and short sequences in the depth of the V1 and V2 pocket; note the abundance of charged residues and N-sites residing in the binding pocket (Fig. 5, Top View). In both structures, the V3 loop, which is imperative for CD4 and coreceptor binding, is buried below the V1–V2 complex (Fig. 5A, B).

As a final point, the greater length variation found in HIV Env sequences from tissues of NPVL subjects may have been present but not sampled in the PVL tissues simply due to the abundance of blood-derived virus. As the NPVL subjects were suppressed until death, it is less likely that these tissues were inundated with HIV contained in recently infected circulating T cells. Furthermore, the SGS approach applied to the NPVL tissues would have likely generated data from all viruses present in the tissues, including archived viruses. We hypothesize that for PVL subject's tissues, whether using either SGS or bPCR, the major blood-associated variants may have overshadowed minor and more diversified variants that were observed in the NPVL subject's tissues.

The increased repertoire of Env structures sampled in NPVL tissues, including both long, short, diversely charged, and glycosylated, fosters intriguing hypotheses concerning persistent tissue-based HIV during cART. Any hypothesis must take into account that tissue-based HIV during cART likely contains archival virus integrated in immune cells, which might be hypermutated or nonfunctional.

One hypothesis suggests that as Env lengths fluctuate over time in blood, the variants repeatedly infect tissue-based immune cells and remain over an extensive period of time with limited evolution, thus the diversity observed in tissue-based HIV is primarily due to its historical nature.

Another hypothesis suggests that tissue-based HIV represents subpopulations of viable and compartmentalized virus that continues to replicate in response to a less understood and complex tissue-based microenvironment. Furthermore, these subpopulations may be masked in subjects who are not on cART due to a massive number of recently replicated viruses.

A third hypothesis combines both ideas, suggesting that the findings presented here are consistent with tissue-based virus populations with both archival and replicating variants. Importantly, functional and longer tissue-based forms of gp120 may be more effective in evading the immune system, as these are tightly linked to immune viral escape and chronic infection^3,14 and thus require attention for consideration as vaccine targets. At the same time, HIV containing functional and shorter Envs in tissues during cART may represent a population of viruses that are more easily transmitted^3,14 or could more easily infect circulating T cells in the absence of cART and thus may play a role in viral rebound post-cART. Thus, NPVL tissue-derived Env sequences, while smaller in numbers than those in PVL tissues, may provide essential mechanistic insights into viral persistence and immune escape that deserves more attention.

As cART has been unable to cure HIV infection, the study of HIV persistence and evolution in tissues is of increased interest. Most studies concerning HIV diversity under cART primarily consider the genes that cART directly targets (e.g., pol, integrase). Remarkably, this study suggests that studies based on Env data derived from tissues of PVL patients may not observe genetic variation representative of an ongoing and dynamic interplay between resident immune cells and tissue-based HIV.

This study also highlights the importance in tissue selection and patient cART history to fully understand tissue-based HIV populations. Using a combination of SGS and deep sequencing approaches to amplify HIV in tissues, such as ultradeep single-molecule real-time sequencing which generates thousands of full-length HIV gene sequences,¹⁷ is extremely useful as the approach can generate long gp120 sequences that reveal the complete sequence Env population in a tissue, including minor variants.