Introduction

Human immunodeficiency virus (HIV-1) rapidly mutates and proliferates in infected individuals. The immune system is a major driver of HIV-1 evolution, as the virus accumulates mutations to escape from host T cells and antibodies ^1–3. Due to the chronic nature of HIV-1 infection, coupled with high rates of mutation and replication, HIV-1 genetic diversity within and between infected individuals is incredibly high. Genetic diversity challenges vaccine development, as vaccine-elicited antibodies must be able to neutralize many strains of the virus to protect against infection ⁴.

However, there exist rare antibodies that are capable of neutralizing a broad range of HIV-1 viruses. These broadly neutralizing antibodies (bnAbs) have therefore been the subject of intense research ^5–8. Eliciting bnAbs through vaccination remains a major goal of HIV-1 vaccine design. However, the development of exceptionally broad antibody responses is rare, and such antibodies typically develop only after several years of infection ^9–12.

Recent years have yielded important insights into the coevolutionary process between HIV-1 and antibodies that sometimes leads to the development of bnAbs. Clinical studies have collected serial samples of HIV-1 sequences from a few individuals who developed bnAbs and characterized the resulting antibodies, their developmental stages, and binding sites ^13–15. The contributions of HIV-1 and its coevolution to bnAb development are complex ^16,17. High viral loads and viral diversity have been positively associated with bnAb development ^16–18. However, superinfection, which can vastly increase HIV-1 diversity, is not always associated with bnAb development ¹⁹, and it does not appear to broaden antibody responses in the absence of other factors ¹⁷.

Here, we sought to characterize the evolutionary dynamics of HIV-1 that accompany the development of bnAbs in clinical data. In particular, we inferred the landscape of selective pressures that shape the evolution of HIV-1 within hosts, reflecting the effects of the immune environment. We first analyzed data from two individuals who developed bnAbs within a few years after HIV-1 infection ^13,14. In both individuals, HIV-1 mutations inferred to be the most beneficial were observed early in infection. In general, mutations that provided resistance to autologous strain-specific antibodies were inferred to be more strongly selected than ones that escaped from bnAbs. We also observed clusters of beneficial mutations along the HIV-1 genome, which were associated with envelope protein (Env) structure.

To confirm the generality of these patterns in a broader sample, we studied recent data from rhesus macaques (RMs) infected with simian-human immunodeficiency viruses (SHIV) that incorporated HIV-1 Env proteins derived from the two individuals above ²⁰. This study also compared patterns of Env evolution in HIV-1 and SHIV in response to host immunity. We observed striking parallels between the inferred fitness effects of Env mutations in RMs and humans, suggesting highly similar selective pressures on the virus despite different host species and differences in individual immune responses. Furthermore, we found that RMs that developed broad, potent antibody responses could clearly be distinguished from those with narrowly focused responses using the evolutionary dynamics of the virus. Specifically, the virus population in individuals who developed greater breadth was distinguished by larger and more rapid gains in fitness than in other individuals. Collectively, these results show high similarity between SHIV evolutionary dynamics in RMs and HIV-1 in humans, and that viral fitness gain is associated with antibody breadth.

Results

Quantifying HIV-1 evolutionary dynamics

We studied HIV-1 evolution accompanying the development of bnAbs in two donors, CH505 and CH848, enrolled in the Center for HIV/AIDS Vaccine Immunology 001 acute infection cohort ²¹. CH505 developed the CD4 binding site-targeting bnAb CH103, which was first detectable 14 weeks after HIV-1 infection ¹³. CH103 maturation was found to be associated with viral escape from another antibody lineage, CH235, that ultimately developed significant breadth ^22,23. CH848 developed a bnAb, DH270, targeting a glycosylated site near the third variable loop (V3) of Env ¹⁴. Similar to the bnAb development process in CH505, escape from “cooperating” DH272 and DH475 lineage antibodies was observed to contribute to the maturation of DH270 (ref. ¹⁴).

To quantify HIV-1 evolutionary dynamics, we sought to infer a fitness model that best explained the changes in the genetic composition of the viral population observed in each individual over time. We assume that the effect on viral fitness of each individual mutation a at each site i is quantified by a selection coefficient s_i(a), with positive coefficients s_i(a) > 0 denoting mutations that are beneficial for the virus and s_i(a) < 0 denoting deleterious ones. We further assume that the cumulative fitness effects of mutations are additive, such that the overall fitness F of a viral sequence is given by the sum of the selection coefficients for all the mutations that it bears. Following this model, we can then quantify the probability of any evolutionary history of the viral population (i.e., the distribution of viral genotypes over time) as a function of the selection coefficients (Methods). We then used MPL ^24,25 to infer the selection coefficients that best fit the viral dynamics observed in data from CH505 and CH848 (Methods). This approach has been successfully applied to study the evolution of HIV-1 (ref. ²⁴), as well as SARS-CoV-2 (ref. ²⁶) and experimental evolution in bacteria ^27,28.

Broad patterns of HIV-1 selection

While the great majority of HIV-1 mutations were inferred to be neutral (s_i(a) ∼0), a few mutations substantially increase viral fitness (Fig. 1). Strongly beneficial mutations occurred in clusters along the genome and preferentially appeared in specific regions of Env (Fig. 2).

Fig. 1.

Fig. 2.

To quantify patterns of selection in Env, we examined the top 2% of mutations inferred to be the most beneficial in CH505 and CH848. The fractions of nonsynonymous mutations within these subsets were 97% and 92% for CH505 and CH848, respectively. These fractions are significantly higher than chance expectations (p = 6.7 ×10⁻³ and 5.6 ×10⁻⁴, Fisher’s exact test; Methods), supporting the model’s ability to accurately infer fitness effects in this data. For CH505, we found 10.9-fold more strongly beneficial mutations in the first variable loop (V1) than expected by chance (p = 2.5× 10⁻³). This is consistent with the presence of V1 mutations conferring resistance to autologous strain-specific antibodies ²². Mutations in V4, a region targeted by CD8⁺ T cells ²², were also 10.0-fold enriched in this subset (p = 9.0 ×10⁻⁵). For CH848, mutations in V1, V3, and V5 were enriched by factors of 14.2, 6.3, and 19.1 among the top 2% most beneficial mutations (p = 5.9 ×10⁻¹⁰, 6.6 ×10⁻³, and 4.8× 10⁻⁶). Mutations in these regions were shown to play a role in resistance to DH270 and DH475 lineage antibodies ¹⁴.

Reversions and mutations affecting N-linked glycosylation motifs were also likely to be beneficial. We define reversions as mutations where the transmitted/founder (TF) amino acid changes to match the subtype consensus sequence at the same site. Among the top 2% most beneficial mutations, reversions were enriched by factors of 19.9 and 17.8 for CH505 and CH848 viruses, respectively (p = 2.1 ×10⁻⁸ and 8.5 ×10⁻¹³), consistent with past work finding strong selection for reversions ^24,30. For CH848, this group also includes several strongly selected mutations observed in Nef (Fig. 1B), a protein that plays multiple roles during HIV-1 infection ^31,32. Mutations affecting N-linked glycosylation motifs (i.e., by adding, removing, or shifting a glycosylation motif) were enriched by factors of 4.6 (p = 7.0 ×10⁻³) and 8. 7(p = 1.4 ×10⁻¹⁰). Changes in glycosylation patterns contributed to antibody escape for both CH505 and CH848 (refs. ^14,22).

Selection for antibody escape

To quantify levels of selection for antibody escape, we computed selection coefficients for mutations that were observed to contribute to resistance to bnAbs as well as bnAb precursors and autologous strain-specific antibodies (Supplementary Tables 1-3) ^22,33–36. We mapped the inferred selection coefficients to the Env protein structure ¹³, highlighting the binding sites for bnAbs and resistance mutations for strain-specific antibodies, as well as important parts of Env (Fig. 2, Supplementary Fig. 1). We also analyzed when these resistance mutations were first observed in each individual.

Overall, we observed stronger selection for escape from autologous strain-specific antibodies and/or changes in glycosylation during the first six months of infection (Supplementary Tables 4-5). This was then followed by more modest selection for escape from bnAb lineages (Fig. 3). For CH505, mutations that conferred resistance to the intermediate-breadth CH235 lineage ²² were less beneficial than top mutations escaping from strain-specific antibodies (Fig. 3A). In turn, resistance mutations for the broader CH103 lineage were less beneficial than CH235 resistance mutations (Fig. 3A). CH848 is similar, with some highly beneficial mutations affecting glycosylation observed early in infection (Fig. 3B). While a few mutations affecting DH270 appear strongly selected, these mutations appeared long before the DH270 lineage was detected (around 3.5 years after infection ¹⁴). Thus, these mutations may have initially been selected for other reasons.

Fig. 3.

Consistent patterns of selection in SHIV evolution

Simian-human immunodeficiency viruses (SHIVs) have numerous applications in HIV/AIDS research ^37,38. Recently, Roark and collaborators studied SHIV-antibody coevolution in rhesus macaques (RMs), which they compared with patterns of HIV-1 evolution ²⁰. Two of the SHIV constructs in this study included envelope sequences derived from CH505 and CH848 transmitted/founder (TF) viruses. There, it was found that 2 out of 10 RMs inoculated with SHIV.CH505 and 2 out of 6 RMs inoculated with SHIV.CH848 developed antibodies with substantial breadth.

To understand whether the patterns of HIV-1 selection observed in CH505 and CH848 are repeatable, and to search for viral factors that distinguish between individuals who develop bnAbs and those who do not, we analyzed SHIV.CH505 and SHIV.CH848 evolution in RMs ²⁰. To prevent spurious inferences, we first omitted data from RMs with <3 sampling times or <4 sequences in total (Methods). After processing, we examined evolutionary data from 7 RMs inoculated with SHIV.CH505 and 6 RMs inoculated with SHIV.CH848 (Supplementary Tables 6-7). We then computed selection coefficients for SHIV mutations within each RM. Reasoning that selective pressures across SHIV.CH505 and SHIV.CH848 viruses are likely to be similar, we also inferred two sets of joint selection coefficients that best describe SHIV evolution in SHIV.CH505- and SHIV.CH848-inoculated RMs, respectively (Methods).

As before, we examined the top 2% of SHIV.CH505 and SHIV.CH848 mutations that we inferred to be the most beneficial for the virus. Overall, we found consistent selection for reversions (17.9- and 14.2-fold enrichment, p = 4.3 ×10⁻¹¹ and 1.2× 10⁻¹¹ for SHIV.CH505 and SHIV.CH848, respectively), with slightly attenuated enrichment in mutations that affect N-linked glycosylation (2.7- and 4.2-fold enrichment, p = 5.4× 10⁻² and 3.4 ×10⁻⁶). However, there is a small subset of mutations that shift glycosylation sites by simultaneously disrupting one N-linked glycosylation motif and completing another, where highly beneficial mutations occur far more often than expected by chance (158.3- and 191.1-fold enrichment, p = 1.7 ×10⁻⁴ and 7.9× 10⁻¹¹ for CH505 and SHIV.CH505, respectively; 118.7-fold and 90.4-fold enrichment, p = 2.0 ×10⁻⁴ and 2.3× 10⁻⁷ for CH848 and SHIV.CH848).

Intuitively, one may expect that strongly beneficial SHIV mutations are more likely to be observed in samples from multiple RMs. However, we found that the number of RMs in which a mutation is observed is only weakly associated with the fitness effect of the mutation (Supplementary Fig. 3). While substantially deleterious SHIV mutations are rarely observed across multiple RMs, neutral and nearly neutral mutations are common. Thus, in this data set, it is not generally true that SHIV mutations observed in multiple hosts must significantly increase viral fitness.

We observed some differences between HIV-1 and SHIV in the precise locations of the most beneficial Env mutations. For example, mutations in V4 are highly enriched in CH505 due to a CD8⁺ T cell epitope in this region, but not in SHIV.CH505 (2.7-fold enrichment, p = 0.13). For SHIV.CH848, beneficial mutations are modestly enriched in V1 and V5 (4.7- and 4.2-fold, p = 4.0 ×10⁻³ and 6.7× 10⁻²), as in CH848, but not for V3 (0.9-fold enrichment, p = 0.22).

Despite some differences in the top mutations, patterns of selection over time in SHIV were very similar to those found for HIV-1. As before, highly beneficial mutations, including ones affecting glycosylation, tended to appear earlier in infection. This was followed by modestly beneficial mutations at later times, including ones involved in resistance to bnAbs in the RMs who developed antibodies with significant breadth (examples in Fig. 4).

Fig. 4.

Detection of SHIV mutations that increase viral load

A major goal of nonhuman primate studies with SHIV is to faithfully recover important aspects of HIV-1 infection in humans. However, due to the divergence of simian immunodeficiency viruses and HIV-1, SHIVs are not always well-adapted to replication in RMs ³⁶. To combat this problem, a recent study identified six SHIV.CH505 mutations that increase viral load (VL) in RMs ³⁶. These mutations result in viral kinetics that better mimic HIV-1 infection in humans.

Our analysis readily identifies the SHIV.CH505 mutations shown to increase VL. Five out of the top six SHIV.CH505 mutations with the largest average selection coefficients are associated with increased VL (Supplementary Table 8). The final mutation identified by Bauer et al., N130D, is ranked tenth. We also find highly beneficial mutations in SHIV.CH848 that are distinct from those in SHIV.CH505 (Supplementary Table 9). Highly-ranked mutations identified here may be good experimental targets for future studies aimed at increasing SHIV.CH848 replication in vivo.

Fitness agreement between HIV-1 and SHIV

Next, we explored the similarity in the overall viral fitness landscapes inferred for HIV-1 and SHIV, beyond just the top mutations. First, we computed the fitness of each SHIV sequence using the joint SHIV.CH505 and SHIV.CH848 selection coefficients inferred from RM data. Then, we computed fitness values for SHIV sequences using selection coefficients inferred from HIV-1 evolution in CH505 and CH848.

We observed a remarkable agreement between SHIV fitness values computed from these two sources (Fig. 5). For both SHIV.CH505 and SHIV.CH848, the correlation between viral fitness estimated using data from humans (i.e., CH505 and CH848) and RMs is strongly and linearly correlated (Pearson’s r = 0.96 and 0.95, p < 10⁻²⁰). This implies that evolutionary pressures on the envelope protein SHIV-infected RMs are highly similar to those on HIV-1 in humans with the same TF Env. In fact, this relationship holds even beyond the SHIV sequences observed during infection. Fitness estimates for sequences with randomly shuffled sets of mutations are also strongly correlated (Supplementary Fig. 4).

Fig. 5.

Evolutionary dynamics forecast antibody breadth

Given the similarity of HIV-1 and SHIV evolution, we sought to identify evolutionary features that distinguish between hosts who develop broad antibody responses and those who do not. Figure 5 shows that SHIV sequences from hosts with bnAbs often reach higher fitness values than those in hosts with only narrow-spectrum antibodies. We hypothesized that stronger selective pressures on the virus might drive viral diversification, stimulating the development of antibody breadth. Past studies have associated higher viral loads with bnAb development and observed viral diversification around the time of bnAb emergence ^13,16,17. Computational studies and experiments have also shown that sequential exposure to diverse antigens can induce cross-reactive antibodies ^39–41.

To further quantify SHIV evolutionary dynamics, we computed the average fitness gain of viral populations in each RM over time. We observed a striking difference in SHIV fitness gains between RMs that developed broad antibody responses and those that did not (Fig. 6). In particular, SHIV fitness increased rapidly before the development of antibody breadth. SHIV fitness gains in RM5695, which developed exceptionally broad and potent antibodies ²⁰, were especially rapid and dramatic. These fitness differences were not attributable to bnAb resistance mutations, which were only moderately selected and generally appeared after bnAbs developed.

Fig. 6.

One outlier in this pattern is RM6072, infected with SHIV.CH505. Antibody development in RM6072 followed a path that was remarkably similar to CH505, including a lineage of antibodies, DH650, directed toward the CD4 binding site of Env ²⁰. However, resistance to the DH650 lineage is conferred by a strongly selected mutation that adds a gly-can at site 234 (T234N, with an inferred selection coefficient of 4.5%). Broadly neutralizing antibodies similar to DH650 are able to accommodate this glycan due to shorter and/or more flexible light chains ^20,42, but DH650 cannot. Antibody evolution in RM6072 thus proceeded along a clear pathway toward bnAb development, but lacked critical mutations to achieve breadth.

Next, we quantified how different types of SHIV mutations contributed to viral fitness gains over time. We examined contributions from VL-enhancing mutations ³⁶, antibody escape mutations ^20,36, other mutations affecting Env glyco-sylation, and reversions to subtype consensus. We found increased fitness gains across all types of mutations in RMs that developed broad antibody responses, compared to those that did not (Supplementary Fig. 5). VL-enhancing mutations, known antibody resistance mutations, and reversions typically made the largest contributions to viral fitness.

Discussion

The induction of bnAbs is a major goal of HIV-1 vaccine design ⁸. Both computational ^39,41,43,44 and experimental ^40,45,46 studies, as well as observations from individuals who developed bnAbs ^13,14,22, suggest that the co-evolution of antibodies and HIV-1 is important to stimulate broad antibody responses. However, a general understanding of the coevolutionary pathways by which bnAbs develop remains incomplete. Here, we quantitatively studied the evolution of HIV-1 and SHIV (featuring HIV-1-derived Env sequences) that accompanies the development of bnAbs.

Our study highlighted how different classes of mutations (e.g., mutations affecting T cell escape or Env glycosylation) affect fitness in vivo. In both HIV-1 and SHIV, we found strong selection for reversions to subtype consensus and some mutations that affected N-linked glycosylation motifs or resistance to autologous strain-specific antibodies. Few CD8⁺ T cell epitopes were identified in this data set, but the T cell escape mutations that we did observe were highly beneficial for the virus. Consistent with past work studying VRC26 escape in CAP256, we observed more modest selection for bnAb resistance mutations ²⁴.

Overall, we found striking similarities between Env evolution in humans and RMs. Importantly, these parallels extend beyond the observation of repetitive mutations, as the number of hosts in which a mutation was observed was only weakly associated with the mutation’s fitness effect (Supplementary Fig. 3). Our inferred Env fitness values in humans and RMs were highly correlated, even for artificial sequences with randomly permuted mutations. This indicates that the evolutionary pressures on Env in HIV-1 and SHIV infection are very similar. Our findings reinforce SHIV as a model system that closely mirrors HIV-1 infection.

We discovered that the speed of SHIV fitness gains was clearly higher in RMs that developed broad antibody responses than in those with narrow-spectrum antibodies. Fitness gains in the viral population preceded the development of bnAbs, and they were not driven by bnAb resistance mutations. This suggests that rapid changes in the viral population are a cause rather than a consequence of antibody breadth. While our sample is limited to 13 RMs and two founder Env sequences, we find a clear separation between RMs that did or did not develop antibody breadth. Thus, the dynamics of viral fitness may serve as a quantitative signal associated with bnAb development.

As noted above, Roark and collaborators also performed a detailed comparison of HIV-1 and SHIV evolution with the same TF Env sequences ²⁰. One of their main conclusions was that most Env mutations were selected for escape from CD8⁺ T cells or antibodies. We found that many antibody resistance mutations identified by Roark et al. are also positively selected in our analysis. Mutations at sites 166 and 169 were shown to confer resistance to a V2 apex bnAb, RHA1, isolated in RM5695 (ref. ²⁰). We inferred moderately positive selection coefficients of 0.49% and 0.43% for R166K and R169K, respectively. The same mutations were found in RM6070, which also developed V2 apex bnAbs, with a selective advantage of 1.7% (Supplementary Table 10). Mutations conferring resistance to autologous strain-specific nAbs were identified at multiple sites by Roark and colleagues: 130, 234, 279, 281, 302, 330, and 334 in RM6072, which developed antibody responses targeting the CD4 binding site (DH650) and V3 (DH647 and DH648) regions. Mutations Y330H and N334S, which confer resistance to V3 autologous nAbs, were detected in all RMs infected with SHIV.CH505, with selective advantages of 3.0% and 4.6% in RM6072, and 1.7% and 3.2% on average across RMs, respectively. Overall, we found that mutations conferring resistance to autologous strain-specific antibodies were common and more strongly selected than bnAb resistance mutations (Supplementary Tables 10-11).

While we observed strong selection for strain-specific antibody resistance mutations, these results could also be affected by the effects of these mutations on viral replication independent of immune escape. In particular, many ssAb resistance mutations are also reversions to the subtype consensus sequence, which have often been observed to improve viral fitness ^24,30. For example, N334S, K302N, and T234N are all reversions. These are among the most beneficial mutations inferred for SHIV.CH505 (Supplementary Table 8). In future work, it would be interesting to attempt to fully separate the fitness effects of mutations due to antibody escape and intrinsic replication (cf. ⁴⁷).

Among other analyses, Roark et al. used LASSIE ⁴⁸ to identify putative sites under selection (Supplementary Tables 12-13). This method works by identifying sites where non-TF alleles reach high frequencies. We found modest overlap between the sites under selection as identified by LASSIE and the mutations that we inferred to be the most strongly selected. For SHIV.CH505, the E640D mutation at site 640 identified by LASSIE is ranked second among 664 mutations in our analysis, and mutations at the remaining 5 sites identified by LASSIE are all within the top 20% of mutations that we infer to be the most beneficial. For SHIV.CH848, the R363Q mutation that is ranked first in our analysis appears at one of the 17 sites identified by LASSIE. Some mutations at the majority of these 17 sites fall within the top 20% most beneficial mutations in our analysis, but some are outliers. In particular, we infer both S291A/P to be somewhat deleterious, with S291P ranked 810th out of 863 mutations.

Beyond the specific context of HIV-1 and bnAb development, our study also provides insight into viral evolution. Despite differences in host biology, the inferred fitness effects of Env mutations in HIV-1 and SHIV were surprisingly similar. This is consistent with the ideas of methods that use sequence statistics across multiple individuals and hosts to predict the fitness effects of mutations ^49–55. However, we found that the relationship between the number of individuals in which a mutation was observed and its inferred fitness effect was fairly weak. This suggests that mutational biases and/or sequence space accessibility may play significant roles in short-term viral evolution, even for highly mutable viruses such as HIV-1 and SHIV. As described above, high frequency mutations were also not necessarily highly beneficial. While the recombination rate of HIV-1 is high, correlations between mutations persist, making it difficult to unambiguously interpret frequency changes as signs of selection ²⁴.

Methods

Data

We retrieved HIV-1 sequences from CH505 (703010505) and CH848 (703010848) from the HIV sequence database at Los Alamos National Laboratory (LANL) ⁵⁶. The rhesus macaque (RM) SHIV sequences ⁵⁷ were obtained from Gen-Bank ⁵⁸. We then co-aligned SHIV.CH505 and SHIV.CH848 sequences with CH505 and CH848 HIV-1 sequences, respectively, using HIValign ⁵⁹.

CH505

CH505 developed two distinct lineages of CD4 binding site (CD4bs) bnAbs, CH103 and CH235 (refs. ^60,61). CH103 antibodies were detectable by 14 weeks after infection and further developed neutralization breadth between 41-92 weeks ⁶². IC50 values of CH235 against the TF virus were 6.5-fold lower than those of CH103 (ref. ⁶⁰). CH235 lineages could neutralize autologous viruses at week 30. However, viruses that acquired mutations at loop D from 53-100 weeks escaped CH235 neutralization ⁶⁰. Although the neutralization breadth of CH235 was not as broad as that of CH103, this lineage played a critical role; escaping mutations from the CH235 lineage stimulated the development of another lineage with broader neutralization depth ⁶⁰. Mutations in loop D enabled the virus to escape from CH235, but these sequential mutations in loop D, such as E275K, N279D, and V281S, favorably bound to the mature CH103 and continuously increased the binding affinity between mature CH103 and loop D ⁶⁰. Gradually, CH103 matured, developing a broader neutralization breadth.

CH848

CH848 developed DH270, a bnAb that targets the glycosylated site adjacent to the third variable loop (V3). DH270 was detectable three and a half years after infection ⁶³. Similar to the CH505 case, the CH848 case exhibited cooperative virus and antibody coevolution. The earlier antibody lineages, DH272 and DH475, could neutralize autologous viruses until week 51 and weeks 15-39, respectively. The virus escaped from DH272 and DH475 afterward, with escape mutations including a longer V1V2 loop. DH270 then developed, with potent and broad neutralization breadth ⁶³.

Rhesus macaques

Chimeric viruses, SHIVs, were constructed by bearing the transmitted/founder (TF) Env from three HIV-1 patients, including CH505 and CH848 (ref. ⁵⁷). In some RMs, SHIV developed similar patterns of mutations to those observed in human donors. In our analysis, we considered RMs with SHIV sequences sampled at at least three points in time. This yielded a set of 7 RMs and 6 RMs for SHIV.CH505 and SHIV.CH848, respectively. The Supplementary Tables 6-7 summarize the number of sequences, time points, and the development of bnAbs for each individual in the SHIV cases as well as HIV-1 cases.

Sequence data processing

Data quality control

To focus our analysis on functional sequences, we removed sequences with more than 200 gaps. To eliminate rare insertions or possible alignment errors, we also masked sites where gaps occurred in more than 95% of sequences within each individual host. To limit errors in virus frequencies, we only considered data from time points with four or more sequences.

Identifying reversions

A mutation is classified as a reversion if the new (mutant) nucleotide matches with the nucleotide at the same site in the HIV-1 consensus sequence from the same subtype. Here, all viruses were subtype C, so we compared with the subtype C consensus sequence as defined by LANL.

Identifying mutations that affect N-linked glycosylation

To identify mutations that affect glycosylation, we search for Env mutations that modify the N-linked glycosylation motif Asn-X-Ser/Thr, where X can be any amino acid except pro-line. We identified three types of mutations affecting glyco-sylation: “shields,” which complete a previously incomplete glycosylation motif, “holes,” which disrupt an existing glyco-sylation motif, and “shifts,” which simultaneously complete one N-linked glycosylation motif and disrupt another.

Enrichment analysis

We used fold enrichment values and Fisher’s exact test to quantify the excess or lack of mutations. For a particular subset of mutations (for example, the top x% beneficial mutations), we first computed the number of mutations in that subset that do (n_sel) and do not (N_sel) have a particular property (for example, nonsynonymous mutations in the CD4 binding site). We then computed the total number of mutations that do and do not have the property (n_null and N_null, respectively) across the entire data set. The fold enrichment value is then . The term quantifies the fraction of mutations having specific properties across the selected mutations, while the denominator is the fraction of all mutations that have the property. Fisher’s exact p values are computed from the 2× 2 table with n_sel, n_null −n_sel in the first row, and N_sel −n_sel, N_null −N_sel − (n_null −n_sel) in the second row ⁶⁴.

Inferring fitness effects of mutations

In this section, we describe the inference framework used to infer the fitness effects of mutations (selection coefficients) from temporal genetic data.

Evolutionary model

We model viral evolution with the Wright-Fisher (WF) model, a fundamental model in population genetics ⁶⁵. In this model, a population of N individuals (viruses or infected cells, in our case) undergo discrete rounds of selection, mutation, and replication. Each individual a has a genetic sequence . Here L is the length of the genetic sequence (number of loci). We write the number of states at each locus (i.e., nucleotides or amino acids) as q. For simplicity, we will use q = 2 and g_i∈ {0, 1} to describe the model below, though we use q = 5 and q = 21 for DNA and amino acid sequences, respectively, in real HIV-1 and SHIV data. In this binary representation, 0 represents a “wild-type” (WT) or TF allele, and 1 is a mutant allele.

We define the fitness of an individual with genetic sequence g by

Here s_i is a selection coefficient, quantifying the fitness effect of a mutant allele at locus i. If s_i > 0, the mutant allele is beneficial (enhancing replication), and if s_i < 0 it is deleterious (impairing replication). By convention, we set the selection coefficient for TF alleles to zero. Individuals with higher fitness values are more likely to replicate than those with lower fitness.

Mutations introduce new genotypes and drive the evolution of the population. Let us define µ^ab as the probability of mutation from genotype a to genotype b per replication cycle. Below, we will express this probability in terms of a mutation rate per site per round of replication. In the analysis of real data, we use asymmetric mutation rates estimated from intra-host HIV-1 data ⁶⁶.

Given these parameters, the WF model describes the dynamics of the frequencies of different genotypes in the population over time. We write the frequency of genotype a at time t as z_a(t). Given that the frequency of genotypes in the population at time t is z(t) = (z₁(t), z₂(t), …, z_M (t)), where M is the total number of genotypes, the probability distribution of the frequency of genotypes in the next generation z(t + 1) is

Here p_a is

where F^a is the fitness value of genotype a, based on Eq. (1). Across K generations, the probability of an entire evolutionary trajectory, defined by the vector of genotype frequencies at each time, is then

Diffusion limit

When the population size is sufficiently large, the evolution of the population defined in Eq. (2) can be reasonably well approximated by a Gaussian process, which is a solution to the Fokker-Planck (forward Kolmogorov) equation ^67,68.

with the drift vector d(t) and the diffusion matrix C(z(t))/N such that

and

Dimensional reduction

While the WF process in genotype space provides valuable insights into genotype dynamics, the mathematical expressions are sometimes challenging to interpret. To obtain more intuitive expressions, we can project the dynamics onto the space of allele frequencies,

One can then find the drift vector d and diffusion matrix C/N in allele frequency space,

and

Here, x_ij is the frequency of individuals with alleles i and j, and is net expected change in frequency of allele i due to mutations, which is given explicitly in Eq. (14) below. The first term in Eq. (9) gives the expected change in the frequency x_i due to the direct fitness effect s_i, while the second term represents the contributions due to indirect or genetic linkage effects with other alleles j.

Maximum path likelihood

Following recent work ⁶⁹, we employed Bayes’ rule to find the selection coefficients that best explain the data. These are the coefficients that maximize the posterior distribution

which is a product of the likelihood of the evolutionary trajectory observed in the data Eq. (4) (under the diffusion limit Eq. (5)) and a prior distribution for the selection coefficients. We chose a Gaussian prior distribution with zero mean and a covariance of I/(Nγ), where I is the identity matrix. This prior distribution penalizes the inferences of large selection coefficients when they are not well-supported by the data. The maximum a posteriori selection coefficients are then given by ⁶⁹

Here C_int, Δx_int, and Δµ_int represent the covariance matrix, vector of frequency changes, and mutational flux integrated over the evolution

The mutational flux µ^flux is characterized by the rates of mutations from nucleotides β to α, denoted by µ_αβ, which are determined from longitudinal HIV-1 populations in untreated patients ⁶⁶. The change of the α nucleotide frequency at locus i due to mutation is given by:

Inverting the integrated covariance matrix effectively reveals the underlying direct allele interactions and resolves the genetic linkage effects.

The shift in the covariance diagonal in Eq. (12), arising from the selection coefficients’ posterior distribution, reflects the uncertainty in the selection distribution. We used γ = 10 for all simulations, but the model is robust to variation in the strength of regularization ⁶⁹. For the mutation rates, we incorporated the transition probabilities among arbitrary DNA nucleotides, estimated from whole-genome deep sequencing of multiple untreated HIV-1 patients followed for 5 to 8 years post-infection ⁶⁶.

Integration of covariance

When the time interval of the observation Δt is sufficiently short, the trajectory of the allele frequency would be continuous and ideally it would be a smooth curve ⁷⁰. To accurately estimate the covariance matrix, we employ piecewise linear interpolation for frequencies. Let τ ∈ [0, 1], then the linear interpolation for a frequency vector can be expressed as:

which yields,

The explicit expression of the integrated covariance is given in refs. ^69,70.

Joint RM model

In addition to fitness models derived from SHIV data for individual RMs, we inferred a joint model under the assumption that virus evolution within each individual RM with the same TF virus is governed by a similar fitness landscape. This method improves inference accuracy from the WF process and deep mutational scanning data ^71–73. The joint path like-lihood for allele frequency trajectories across RMs with the same TF virus is then

Here, x^α is the allele frequency of the α-th individual and R is the number of replicate individuals (i.e., the number of RMs sharing the same TF virus). The initial state is p(x^α(t₀)) = p(x(t₀)) for all α for individuals with the same TF virus. The solution of the joint path likelihood is given by

Here, the overbar denotes the sum over the replicate RMs. We emphasize that the joint selection coefficients in

Eq. (17) are not the same as selection coefficients that are simply averaged across RMs with the same TF virus. The joint selection coefficients are more robust, as they are guided by the level of evidence within each individual rather than naive averaging.

Geometrical interpretation of the fitness comparison

The Pearson values we utilized to compare the fitness landscapes denoted as F_s(g) := F (g |s) and F_h(g) := F (g |h) can be expressed by the following simple relation:

Here, Cov(F_s, F_h) and Var(F_s) represent the covariance and variance values estimated from the samples being compared, (F_s(gⁿ), F_h(gⁿ))_n. C is the covariance matrix defined between arbitrary loci. The last equation can be interpreted as an angle between two vectors, s and h, with a metric matrix C; if s = h, the Pearson value clearly becomes 1. However, the “similarity” also depends on how these vectors are projected by C; eigenmodes associated with larger variance of statistics will be more emphasized. The last expression readily implies an interpretation for the case of shuffled sequences; shuffling the sequences equates to diluting the covariance between loci, resulting in the metric matrix becoming a diagonal matrix. Removing the off-diagonal elements corresponds to lifting the constraints on the fitness landscape.

Acknowledgements

The work of K.S.S. and J.P.B. reported in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award Number R35GM138233.

Additional information

Author contributions

All authors contributed to the research design, methods development, interpretation of results, and writing the paper. K.S.S. performed simulations and computational analyses. J.P.B. supervised the project.

Code availability

Sets of processed data and computer code used in this study are available in the GitHub repository of https://github.com/bartonlab/paper-HIV-coevolution. This repository contains source files that process HIV-1 and SHIV sequences, infer selection coefficients, and identify and characterize mutations. The included Jupyter note-books can be run to reproduce the figures presented here. The original HIV-1 sequences can be retrieved from the LANL database (https://www.hiv.lanl.gov/content/index), and SHIV sequences can be found at GenBank (https://www.ncbi.nlm.nih.gov/).