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 sitetargeting 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. In recent years, a wide variety of approaches have been developed to infer the fitness effects of mutations from temporal genetic data ^24–43. The vast majority of these methods focus on a single locus at a time, ignoring correlations between genetic variants at different loci. While the recombination rate of HIV-1 is high ^44,45, the virus also evolves under strong selection, which can lead to interference between clones with different beneficial mutations ^39,46–50. Thus, we applied MPL, an inference method that systematically accounts for genetic correlations ^39–43, to estimate the fitness effects of HIV-1 mutations.

Model overview

Here, we provide a brief overview of the key steps in the MPL approach to inferring selection. Further details are available in Methods and in prior work ^39–43. First, 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. That is,

where represents the viral sequence, with qual to one if genotype α has allele a at site I and zero otherwise. L is the length of the genetic sequence, and q represents the number of genetic states (i.e., q = 5 for nucleotide sequences and 21 for amino acids, including gaps/deletions).

We assume that viral replication is stochastic, where viruses with higher fitness are more likely to spread infection to new cells than ones with lower fitness. Let us write the number of viruses of each genotype in the population at time t as n(t) = (n₁(t), n₂(t), …, n_M (t)), where is the total number of possible genotypes. In our model, the probability of obtaining a new distribution of genotypes n(t + 1) in the next generation is multinomial,

with N = ∑_αn_α the total population size. The probabilities p_α(n(t)) are influenced by fitness, as well as mutation, recombination, and the current frequency of each genotype (see Methods). Essentially, viruses that are fitter than the current population average are more likely to increase in frequency, while ones that are less fit are likely to decline. Mutations and recombination introduce genetic variation into the population. The population size N determines the relative stochasticity of the dynamics, with smaller populations having more random fluctuations than larger ones.

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). Assuming that the population size N is large and the selection coefficients are small (such that |s| ≪ 1), we can write an analytical expression for the selection coefficients that best fit the viral dynamics observed in data ^39,40:

Here C_int is the covariance matrix of allele frequencies (i.e., the linkage disequilibrium matrix) integrated over time, and γ is a regularization parameter. The terms Δx_int and Δu_int represent the net change in allele frequency and the total expected change in allele frequency due to spontaneous mutations alone, respectively (see Methods for details). Intuitively, this expression says that net allele frequency changes that cannot be explained by mutations are likely due to selection, either on that specific allele or the associated genetic background, which is quantified by C_int. Alleles that have large, rapid changes in frequency are more likely to be under strong selection than those with smaller, slower frequency changes. Beyond HIV-1 (ref. ³⁹), this approach has been successfully applied to study the evolution of SARS-CoV-2 (ref. ⁴³) and experimental evolution in bacteria ^51,52.

Broad patterns of HIV-1 selection

As described above, we used MPL ^39,40 to infer the selection coefficients that best fit the viral dynamics observed in data from CH505 and CH848. 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 ¹⁴. To test whether our results might be biased by overall sequence variability, we examined the relationship between our inferred selection coefficients and entropy, a common measure of sequence variability. Overall, we found only a modest correlation between selection and entropy, suggesting that the signs of selection that we observe are not due to increased sequence variability alone (Supplementary Fig. 1).

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 ^39,54. 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 ^55,56. 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,57–60. 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. 2). 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 ^61,62. 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. 4).

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 immun-odeficiency 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. 5). In contrast, the inferred fitness landscapes of CH505 and CH848, which share few mutations in common, are poorly correlated (Supplementary Fig. 6). This suggests that the similarities between viral fitness values in humans and RMs are not artifacts of the model, but rather stem from similarities in underlying evolutionary drivers.

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 ^63–65.

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 glycan 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,66, 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,60, other mutations affecting Env glycosylation, 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. 7). VL-enhancing mutations, known antibody resistance mutations, and reversions typically made the largest contributions to viral fitness.

Robustness of inferred selection to changes in the fitness model and finite sampling

In the analysis above, we used a simple model where the net fitness effect of multiple mutations is simply equal to the sum of their individual effects. Recently, methods have also been developed that can infer epistatic fitness effects from data, which include pairwise interactions between mutations ^40,41. We reanalyzed these data to examine how inferred fitness changes when epistasis is included in the model, using the approach of Shimagaki and Barton ⁴¹ (Methods). Overall, the inferred epistatic interactions were modest (Supplementary Fig. 8). In CH505, we found that the CD4 binding site, V1 (especially sites 136–146 in HXB2 numbering) and V5 regions were modestly but significantly enriched in the most beneficial (top 1%) of epistatic interactions (2.5-, 1.2- and 1.8-fold enrichment with p = 1.0 × 10⁻²¹, 6.3 × 10⁻⁶ and 6.3 × 10⁻⁵, respectively). Epistatic interactions between N280S/V281A and E275K/V281G, which confer resistance to CH235 (ref. ²²), ranked in the top 6.5% and 13.0% of interactions. In CH848, we found 1.3-, 1.5-, and 2.3-fold enrichment in strong beneficial epistatic interactions in the CD4 binding site, V4, and V5 regions, respectively (p = 4.0 × 10⁻⁶, p = 2.5 × 10⁻¹⁴ and 3.2 × 10⁻¹⁹).

To compare the typical fitness effects of individual mutations in the model with epistasis to those in the additive model, we computed effective selection coefficients for the epistatic model. For each mutant allele a at each site i, we computed the average difference in fitness between sequences in the data set with the mutation and hypothetical sequences that are the same as those in the data, except with the mutant allele a reverted to the TF one. In this way, the effective selection coefficient measures the typical effect of each mutation in the data set, while also accounting for epistatic interactions with the sequence background. We found that the effective selection coefficients were highly correlated with the selection coefficients from the additive model (Supplementary Fig. 9). We also found strong agreement between the additive and epistatic model fitness values for each sequence in both HIV-1 and SHIV data (Supplementary Fig. 10).

Finite sampling of sequence data could also affect our analyses. To further test the robustness of our results, we inferred selection coefficients using bootstrap resampling, where we resample sequences from the original ensemble, maintaining the same number of sequences for each time point and subject. The selection coefficients from the bootstrap samples are consistent with the original data (see Supplementary Fig. 11 for a typical example), with Pearson’s r values of around 0.85 for HIV-1 data sets and 0.95 for SHIV data sets, respectively.

Discussion

HIV-1 evolves under complex selective pressures within individual hosts, balancing replicative efficiency with immune evasion. Here, we quantitatively studied the evolution of HIV-1 and SHIV (featuring HIV-1-derived Env sequences) across multiple hosts, including some who developed broad antibody responses against the virus. 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: the number of hosts in which a mutation was observed was only weakly associated with the mutation’s fitness effect (Supplementary Fig. 4). Our inferred Env fitness values in humans and RMs were highly correlated, indicating that the functional and immune constraints shaping Env evolution in HIV-1 and SHIV infection are very similar. Our findings therefore 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.

The induction of bnAbs is a major goal of HIV-1 vaccine design ⁸. Both computational ^63,65,67,68 and experimental ^64,69,70 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. Our results could thus inform HIV-1 vaccine research. While precise immune responses and viral escape pathways can differ across individuals, the quantitative similarity in viral evolutionary constraints across humans and RMs suggests that SHIV data can provide a valuable source of information about Env variants that contribute to bnAb development, especially when detailed longitudinal data from humans does not exist. While the concept of sequential immunization is well-established ^{10,63,64,71,72}, our findings also suggest a possible new design principle. Immunogens could be engineered to reproduce the dynamics of viral population change that are associated with rapid fitness gains, which we found to precede the emergence of bnAbs. This emphasis on broader, population-level dynamics could complement investigations of the molecular details of virus and antibody coevolution.

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).

We note that our conclusions about the phenotypic effects of HIV-1 mutations under selection are constrained by the available data. 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 ^39,54. 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 anti-body escape and intrinsic replication (cf. ⁴²). Although we have systematically compiled information about mutations known to affect antibody resistance and glycosylation, this data is necessarily incomplete. Some of the strongly beneficial mutations with unknown functional effects that we observe could therefore reflect escape from unmapped immune responses.

There are additional methodological and technical limitations that should be considered in the interpretation of our results. Most notably, we assume that the viral fitness landscape is static in time. While we do not expect selection for effective replication (“intrinsic” fitness) to change substantially over time, pressure for immune escape could vary along with the immune responses that drive them. In prior work, we have found that constant selection coefficients typically reflect the average fitness effect of a mutation when its true contribution to fitness is time-varying ^42,43. This may not adequately description mutational effects that undergo large or rapid shifts in time. Future work should also examine temporal patterns in selection for individual mutations.

While we found a strong relationship between viral fitness dynamics and the emergence of bnAbs, it may not be true that the former stimulates the latter. For example, bnAbs may have been present within each host before they were experimentally detected. Rapid viral fitness gains within hosts that developed broad antibody responses could then have been driven by undetected bnAb lineages. However, we did not find strong selection for known bnAb resistance mutations, and in at least one case (RM5695), rapid fitness gains (roughly 2 weeks after infection) substantially preceded bnAb detection (16 weeks). Still, given the limited size of the data set that we studied, it is unclear the extent to which our results will transfer to larger and broader data sets.

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 across hosts and related host species. Parallels between the HIV-1 and SHIV fitness landscapes that we infer suggest that there are strong constraints on viral protein function, with few paths to significantly higher fitness. This is consistent with the ideas of methods that use sequence statistics across multiple individuals and hosts to predict the fitness effects of mutations ^74–80. However, 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 ³⁹.

Our results also point to strong similarities in the immune environment across closely related host species, including preferential targeting of specific parts of viral surface proteins by antibodies. This is supported by the enrichment of beneficial mutations within variable loop regions and at sites that affect the glycosylation of Env. However, despite these constraints, there may still exist a large number of neutral or nearly-neutral mutational paths that remain unexplored.

Overall, our findings support the potential predictability of viral evolution, at least over short time scales. While there are contingencies in evolution – for example, disparate host immune responses or strong epistatic constraints between mutations – these are not so pervasive that they completely change the effective viral fitness landscape or paths of evolution across hosts, given the same founder virus sequence. Similar observations of parallel evolution in HIV-1 have been reported in monozygotic twins infected by the same founder virus ⁸¹, common patterns of immune escape across hosts ^78,82 and drug resistance ^83–85, and long-term experimental evolution ⁸⁶. Our results thus contribute to a growing body of research identifying predictable features in viral evolution. Understanding such features could ultimately inform practical applications such as anticipating the emergence of drug resistance or designing vaccines to limit likely pathways of escape.

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. ^91,92). 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 proline. We identified three types of mutations affecting glycosylation: “shields,” which complete a previously incomplete glycosylation motif, “holes,” which disrupt an existing glycosylation 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 genotype α is represented by a sequence , where is equal to one if genotype α has allele a at locus i and zero otherwise. Here L and q represent the length of the genetic sequence (number of loci) and the number of statues at each locus (i.e., number of nucleotides or amino acids), respectively. We use q = 5 and q = 21 for DNA and amino acid sequences, respectively, in real HIV-1 and SHIV data.

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

Here s_i(a) is a selection coefficient, quantifying the fitness effect of allele a at locus i. If s_i(a) > 0, the allele a is beneficial (enhancing replication), and if s_i(a) < 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 µ^αβ as the probability of mutation from genotype α to genotype β 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 α at time t as z_α(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_α is

where F^α is the fitness value of genotype α, based on Eq. (4). 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. (5) can be reasonably well approximated by a Gaussian process, which is a solution to the Fokker-Planck (forward Kolmogorov) equation ^98,99.

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(a, b) is the frequency of individuals with alleles a and b at loci i and j, and u_i(a) is net expected change in frequency of allele a at i due to mutations, which is given explicitly in Eq. (17) below. The first term in Eq. (??) gives the expected change in the frequency x_i(a) due to the direct fitness effect s_i(a), 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. (7) (under the diffusion limit Eq. (8)) 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 Δu_int represent the covariance matrix, vector of frequency changes, and mutational flux integrated over the evolution

The mutational flux u is characterized by the rates of mutations from nucleotides b to a, denoted by µ_ab, which are determined from longitudinal HIV-1 populations in untreated patients ⁹⁷. The change of the a 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. (15), arising from the selection coefficients’ posterior distribution, reflects the uncertainty in the selection distribution. We used γ = 10 for all data sets, 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,

For simplicity in notation, we omitted nucleotide indices. The explicit expression of the integrated covariance is given in refs. ^100,101.

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 ^102–104. The joint path likelihood for allele frequency trajectories across RMs with the same TF virus is then

Here, x^r is the allele frequency of the r-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^r(t₀)) = p(x(t₀)) for all r 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. (20) 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.

Epistatic fitness model

We consider the following pairwise epistatic fitness function, which depends on epistatic interactions s_ij across all possible pairs of loci:

Our goal is to obtain the epistatic interactions s_ij(a, b) as well as the selection coefficients s_i(a) from temporal genetic sequences.

The basic logic for inferring these fitness parameters parallels the additive case. The only practical difference is that epistatic interactions can influence the dynamics of additive and pairwise mutation frequencies. Under the diffusion limit, we obtain the following drift terms, which align with Eq. (12) in the additive model,

and diffusion matrices,

Here, u_i and u_ij represent the expected frequency changes due to mutations for additive and pairwise terms, while v_i represents the changes in pairwise frequencies due to recombination. These explicit expressions indicate that the u_i remains the same as in the additive fitness case; therefore, Eq. (17) holds. The pairwise term is given as

and the v_ij is expressed as

where r denotes the recombination rate. In this study, we set r = 10⁻⁵. More detailed derivations are provided in the previous studies ^102,103.

The technical challenge of the epistasis inference is that the diffusion matrix C involves third- and fourth-order interactions, and the number of matrix elements scales as 𝒪 ((qL)⁴), while the computational cost to invert it scales as 𝒪 ((qL)⁶). Recently, a more efficient computational method was proposed, reducing both the necessary memory usage and computational times by 𝒪 ((qL)²) (ref. ¹⁰³). The essential idea involves factorizing the higher-order covariance matrix using the rectangular matrix Ξ∈ ℝ^D×d such that C = ΞΞ^⊤, where D scales as 𝒪 ((qL)²) while d ≪ D. This method resolves the linear equation without obtaining an explicit expression of C and avoids any computations involving more than O((qL)²) operations ¹⁰³.

Gauge transformation

Since constant shifts in fitness parameters do not affect relative fitness, it is always possible to transform fitness values without changing the resulting genotype distribution. For example, in the additive model, individual selection coefficients can be shifted as

where is the allele of a chosen reference genotype (e.g., the TF sequence) at site i. This transformation preserves relative fitness but it ensures that , making inferred selection coefficients more interpretable.

In statistical physics, such transformations, where model parameters are changed without altering the underlying probability distribution—are referred to as gauge transformations ^105,106. Similar transformations have been employed in recent studies to improve interpretability and sparsity in epistatic models ¹⁰⁷. We can apply an analogous transformation to epistatic interactions:

Under this transformation, any selection coefficients or epistatic terms involving TF alleles are zero by definition, while relative fitness remains unchanged.

Regularization

Regularization is used to reduce the effective number of parameters in the fitness model. In our analysis, we applied strong regularization (γ = 10¹⁰) to any selection or epistatic coefficients involving TF alleles, ensuring they are effectively zero under the gauge transformation. Following prior work ¹⁰⁸, we penalized epistatic interactions between loci more than 50 nucleotides apart on the reference sequence with the same strong regularization. We also used the same moderate regularization value of γ = 50 for all other epistatic terms, and used γ = 10 for selection coefficients, consistent with the additive model.

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/).

Funding

National Institutes of Health (R35GM138233)