The replication-competent HIV reservoir is a genetically restricted, younger subset of the overall pool of HIV proviruses persisting during therapy, which is highly genetically stable over time

We reconstructed the within-host HIV evolutionary histories of seven women who initiated ART at a median of 9 (range 1.9–12) years following their estimated infection dates, which were calculated as the midpoint between their last negative and first seropositive HIV visits (Fig. 1; Table 1). We leveraged this information to characterize the genetics and dynamics of proviruses sampled longitudinally during ART, as well as HIV RNA sequences rebounding in plasma. Together, we analyzed 1,092 single-genome-amplified intact HIV RNA env-gp120 sequences (median 181, range 50–239 per participant) from a median of 9 (range 2–13) plasma samples collected over a median of 8.3 (range 0.8–11.8) years pre-ART, where these sequences were previously published for three participants (32) (Fig. 1). In addition, we analyzed 926 intact proviral env-gp120 sequences (median 150, range 42–182 per participant) from a median of 3 time points (range 1–4) spanning a median of 8.7 (range 2.8–12.3) years during ART. For participant 1, we also analyzed 114 plasma HIV RNA env-gp120 sequences isolated during a viremia rebound event, and for participant 5, we also analyzed four plasma HIV RNA env-gp120 sequences isolated during initial loss of viral control. All participants had HIV subtype B, with no evidence of dual or super-infection (Fig. S1). As expected (8, 32–34), the overall extent of within-host HIV diversity correlated strongly with the duration of untreated infection (Spearman’s ρ = 0.85, P = 0.03; Fig. S1, inset).

For participants 1 through 6, we sampled proviruses at a minimum of three time points during the initial years of ART, allowing us to investigate clonal dynamics. The overall percentage of putatively clonal sequences, defined as those that matched at least one other sequence with 100% nucleotide identity in env-gp120, ranged from 9% (participant 4) to 58% (participant 3) (Fig. 2A). In all participants, we recovered clones that were observed at only one on-ART time point (gray slices, Fig. 2B) and across more than one time point (colored slices, Fig. 2B). In participants 3, 5, and 6, we recovered at least one clone that persisted across all on-ART time points. The proportion of clonal proviruses increased or remained stable over time in five of six participants (Fig. 2C).

We used a phylogenetic approach (7) to investigate the ages, within-host evolutionary origins, and dynamics of proviral lineages persisting during ART. Only intact, non-hypermutated sequences that showed no evidence of within-host recombination were included in this analysis (recombinants were analyzed separately; see below). To mitigate the inherent uncertainty in within-host HIV evolutionary reconstruction and to allow us to estimate error in the parameters of interest (e.g., proviral integration dates; population genetic structure), we inferred a minimum 1,500 trees per participant and conditioned results across all trees. We rooted each tree at the location that maximized the correlation between the root-to-tip distances of the pre-ART plasma HIV RNA sequences and their sampling dates (as within-host sequence divergence from the transmitted/founder virus increases over time during untreated infection (32, 34, 35)). This root location represents the most recent common ancestor of the data set, which should be the transmitted/founder virus (or a close descendant) in this cohort of seroconverters. We then fit a linear model to each tree relating the root-to-tip genetic distances of distinct pre-ART plasma env-gp120 sequences to their sampling dates. Here, the slope represents the within-host pre-ART env-gp120 evolutionary rate, and the x-intercept represents the phylogenetically estimated infection date. This linear model was then used to convert the root-to-tip distance of each post-ART sequence of interest to its integration date. Each sequence’s estimated integration date was then averaged over all within-host trees that passed quality control (QC) (see Materials and Methods) and reported along with its 95% highest posterior density (HPD) interval.

Participant 1 was estimated to have acquired HIV in late 1995 but only initiated ART in January 2008 (Fig. 1 and 3A). We inferred 15,000 within-host phylogenies relating 207 plasma HIV RNA env-gp120 sequences collected over 12 years of untreated infection (these sequences were originally published in (32)), 97 plasma HIV RNA env-gp120 sequences collected at three time points during a viremia rebound event, and 90 proviral env-gp120 sequences collected at four subsequent time points after viremia was re-suppressed on ART (Table 1; Table S1). Of these trees, 7,218 passed QC (see Materials and Methods and Table S1). The trees exhibited the ladder-like shape typical of within-host HIV evolution, where plasma HIV sequences sampled during untreated infection were increasingly divergent from the root over time (example phylogeny in Fig. 3B). This shape is the result of serial genetic bottlenecks imposed by immune responses from which the virus continually escapes; the selective sweeps that characterize this process can be seen in the adjacent amino acid highlighter plot. The linear model relating the root-to-tip distances of pre-ART plasma HIV RNA sequences in this tree to their collection dates is shown in Fig. 3C, where the slope of this line represents the within-host pre-ART env-gp120 evolutionary rate and the x-intercept represents the phylogenetically estimated infection date. Conditioning these metrics over all QC-passed trees revealed a mean root date of December 1995 (95% HPD, from August 1995 to March 1996), consistent with the participant’s clinically estimated infection date, and a mean env-gp120 evolutionary rate of 4.8 × 10⁻⁵ (95% HPD 3.8 × 10⁻⁵ – 6.0 × 10⁻⁵) substitutions per nucleotide site per day (Table S1). All subsequent analyses were restricted to distinct HIV sequences per time point, as our goal was to examine the diversity and ages of distinct proviral sequences in a way that is not influenced by clonal expansion.

Beginning in October 2010, participant 1 experienced a prolonged viremia event that occurred after a documented change in ART regimen. Unfortunately, neither adherence nor drug resistance data were available to confirm the cause of the rebound. Nevertheless, this event allowed us to investigate the within-host evolutionary origins of HIV sequences rebounding in plasma, which we sampled at three time points between October 2010 and October 2011. All sequences sampled at the initial time point interspersed with plasma sequences that circulated in the 3 years prior to ART (Fig. 3B), where our phylogenetic methods dated these sequences to between late 2005 and the time of sampling (Fig. 3C and D). Rebound sequences subsequently sampled in April and October 2011 also dated to this era, which is expected given that many of these sequences were descendants of the initial rebounding population (Fig. S2 shows an enlarged portion of the phylogeny, where likely descendants are indicated by green arrows). Nevertheless, a monophyletic clade of closely related rebound sequences isolated in April 2011 was less divergent from the root and distinct from the rebound sequences sampled the previous October (Fig. S2, black bracket). This observation was consistent across all trees (see consistently earlier integration dates for these sequences in Fig. 3D). This suggests that these were not descendants of the initial rebounding population but rather descendants of a provirus or clonal infected cell population that reactivated independently near this time.

In striking contrast to the rebound HIV sequences, participant 1’s distinct proviral sequences subsequently sampled on ART in 2013, 2015, 2017, and 2018 interspersed throughout the entire phylogeny, nearly all the way back to the root (Fig. 3B). Though the inferred integration dates of these proviruses were slightly skewed toward the years leading up to ART, as well as the subsequent ART interruption period, numerous older proviruses were recovered, including one dating to September 1996, less than a year following infection (Fig. 3D). The integration date distributions of distinct proviruses sampled longitudinally on ART were highly stable over time (Kruskal-Wallis P = 0.6; Fig. 3D). Notably, these proviruses were overall significantly older than those that had previously rebounded in plasma (Mann-Whitney P = 0.0005 when comparing the ages of all proviruses vs. those of the initially rebounding population; P < 0.0001 when comparing all proviruses to the entire rebounding population). We also observed proviruses that were identical or near identical to sequences that had previously rebounded in plasma (Fig. S2, black arrows), consistent with reservoir re-seeding during this rebound event.

We next investigated whether the distinct proviral sequences sampled during ART showed any evidence of changing population genetic structure over time. This could happen, for example, if distinct within-host lineages were being eliminated during the initial years of ART or if proviral populations sampled at different time points were distinct from one another (which could occur, for example, if reservoir cells re-entered blood from a compartmentalized tissue population). Neither test revealed significant evidence of population structure: AMOVA yielded P = 0.3, while the CC test, conditioned over all passing trees, yielded a mean P = 0.16 (Fig. 3E). Therefore, even though proviral clonality increased during the first 9 years of ART (Fig. 2C), the number and composition of distinct proviral lineages remained consistent over time. In contrast, plasma rebound HIV sequences were significantly compartmentalized compared with the overall proviral pool (P = 0 for both AMOVA and CC), indicating that only a restricted subset of persisting proviruses re-seeded viremia at that time.

Within-host recombinants, which are an important source of within-host HIV diversity (36, 37), cannot be dated using the phylogenetic approaches used here. For recombinants, we instead used RDP4 (38) to obtain the origin date of each parental sequence component. We were particularly interested in the possibility of recombination between sequences from different infection eras, a phenomenon that has been predicted by mathematical modeling (39). The existence of such recombinants would support the notion that reservoir cells can be superinfected with HIV from another within-host era.

In participant 1, we identified 40 distinct recombinant env-gp120 proviral sequences that were sampled on ART between 2013 and 2018 (Fig. 4). None contained a parent sequence from the first 7 years of infection. Rather, 5 were mosaics of sequences that rebounded during the 2010–2011 treatment interruption, 6 were mosaics of sequences that circulated in the 3 years prior to ART, and 22 were mosaics of rebound sequences and pre-ART sequences from 2005 to 2007, which fueled the rebound event. The final five were sequences whose minor parent could not be identified. Two mosaics comprising 2003 and rebound sequences were observed, but these parent sequences could plausibly have circulated at the same time, as the root-to-tip divergences of sequences from these periods overlapped one another (see Fig. 3C). As such, there was no convincing evidence in this participant of recombinants between viruses from substantially different infection eras. Similarly, the 17 recombinants identified among the 2010–2011 plasma rebound sequences were all mosaics of rebound sequences and/or pre-ART sequences that fueled the rebound (Fig. 4).

Participant 2 was estimated to have acquired HIV in January 2003 (Table 1). Though ART was briefly initiated in 2007, durable suppression was not achieved until ART was re-initiated in January 2012 (Table 1; Fig. 5A). We inferred 5,250 phylogenies from 239 plasma HIV RNA env-gp120 sequences isolated during untreated infection, along with 61 proviral env-gp120 sequences sampled at 2, 6, and 7 years post-ART. Of these, 3,842 trees passed QC, yielding a mean root date of May 2003 (95% HPD, from January to September 2003; Table S1). Proviral sequences sampled during ART interspersed throughout nearly the whole tree (example phylogeny and root-to-tip divergence plot in Fig. 5B and C), where the oldest provirus dated to late 2004, approximately 2 years after infection (Fig. 5D). The integration date distributions of distinct proviruses sampled over 7 years on ART remained consistent over time (Kruskal-Wallis P = 0.3; Fig. 5D) and showed no evidence of temporal population structure (AMOVA P = 0.2; CC mean P = 0.1; Fig. 5E).

Four recombinant proviral sequences were recovered from participant 2, three of which were mosaics of sequences that circulated in the 2 years prior to ART (Fig. S3A). In contrast, one recombinant had a 5′ half that dated to 2006 and a 3′ half that dated to 2010 (Fig. S3A; asterisk). Of note, the root-to-tip divergence measurements of the plasma HIV RNA sequences isolated in 2006 and 2010 did not overlap one another (see Fig. 5C), suggesting that sequences with such distinct divergence measurements could not have co-circulated during ongoing HIV evolution. Instead, the recovery of this recombinant sequence suggests that a reservoir cell dating to 2006 became superinfected with a sequence circulating in 2010 (or that a cell infected in 2010 became superinfected with HIV reactivated from 2006). This superinfected cell would then have needed to produce an infectious virion carrying both parent genomes, which then yielded a novel recombinant provirus upon infection of a new cell.

Participant 3 was estimated to have acquired HIV in July 2002 and initiated ART in January 2008 (Table 1; Fig. 6A). We inferred 4,500 phylogenies from 140 pre-ART plasma env-gp120 sequences, published in (32), along with 96 proviral sequences sampled after 4, 7, and 9 years on ART. All passed QC, yielding a mean root date of November 2001 (95% HPD, July 2001 to March 2002), which was slightly earlier than the clinically estimated infection date of July 2002 (Table S1). Again, distinct proviral sequences sampled on ART interspersed throughout the tree (example phylogeny and divergence plot shown in Fig. 6B and C). We also observed a number of expanded clones. These included one relatively near the root that we recovered at both the second and third on-ART sampling time points (Fig. 6B, black arrow) and another more divergent one that we recovered at all on-ART time points (Fig. 6B, green arrow). The oldest recovered provirus, isolated in 2014, was estimated to have integrated in May 2002 (Fig. 6D). The integration date distributions of distinct proviruses were stable during ART (Kruskal-Wallis P = 0.9; Fig. 6D), with no evidence of changing population structure (AMOVA P = 0.1; CC mean P = 0.25; Fig. 6E). All 12 recombinant proviral sequences collected during ART were mosaics of sequences that circulated in the 3 years prior to ART (Fig. S3B).

Participant 4 was estimated to have acquired HIV in July 1995 and initiated suppressive ART in June 2006 (Table 1; Fig. 7A). We inferred 6,000 within-host phylogenies from 195 pre-ART plasma env-gp120 sequences along with 141 proviral env-gp120 sequences sampled at four time points up to 12 years post-ART (Table 1). All trees passed QC (example phylogeny and divergence plot in Fig. 7B and C), yielding a mean estimated root date of February 1995 (95% HPD, from October 1994 to June 1995; Table S1), which was only slightly earlier than the clinically estimated infection date. Proviruses isolated during ART had integrated throughout the entire course of untreated infection and included five distinct sequences, isolated in March and September 2018, that dated to 1995 (i.e., the first 6 months after transmission). In contrast to participants 1–3, the age distributions of distinct proviruses sampled in March and September 2018 (~12 years after ART initiation, the longest follow-up of any participant in the study) were on average older than those sampled in the earlier years of ART (Kruskal-Wallis P = 0.0004; pairwise post-test P-values shown in Fig. 7D). Both AMOVA (P = 0.001) and CC (mean P = 0.004) also supported a change in proviral population structure over this period (Fig. 7E).

Eleven recombinant proviral sequences were recovered from participant 4, 10 of which were mosaics of sequences that circulated in the 3 years prior to ART initiation (Fig. S4A). Of note, and similar to participant 2, one recombinant was a mosaic of a sequence from the year 2000 and another from 2005 (Fig. S4A; asterisk). The root-to-tip divergence measurements of plasma HIV RNA sequences circulating in 2000 and 2005 would not have been expected to overlap (see Fig. 7C), again supporting the notion that this recombinant arose as a result of a reservoir cell becoming superinfected with HIV from different infection eras.

Participant 5 acquired HIV around March 2008 and initiated ART in February 2010 (Table 1; Fig. 8A). We inferred 1,500 phylogenies relating 50 env-gp120 sequences collected at two pre-ART time points, 48 proviral env-gp120 sequences collected at 2 and 5 years after ART initiation, 4 HIV RNA env-gp120 sequences collected in year 8 when viremia initially rose to 295 HIV RNA copies/mL, and 26 proviral sequences collected in year 9, after viral control was lost (Table 1). Unfortunately, neither drug resistance nor adherence data were available to confirm the cause of the viremia control loss, and there were no documented ART regimen changes around this time. All phylogenies passed QC (Table S1) and yielded a mean estimated root date of December 2008 (95% HPD, October 2008 to February 2009), suggesting that we did not reconstruct all the way back to transmission (example phylogeny and divergence plot in Fig. 8B and C). Proviruses sampled during ART integrated throughout untreated infection, with the earliest dating to May 2009, about a year after transmission. In contrast, and similar to the observations from participant 1, both distinct plasma HIV RNA sequences isolated during the initial loss of viremia control dated to around ART initiation (Fig. 8B). Although the age difference between the emerging plasma viruses and the overall proviral pool did not reach statistical significance due to small numbers (P = 0.11), the emerging plasma viruses were clearly younger (Fig. 8D). The integration date distributions of proviruses sampled after 2, 5, and 9 years of ART were not significantly different from one another (Kruskal-Wallis P = 0.1), though the year 9 proviral pool featured a slightly larger number of younger proviruses than the first two time points, likely due to the loss of viremia control at this time (Fig. 8D). The P-values of the population structure tests suggested a similar shift (AMOVA yielded P = 0.02 and CC yielded a mean P = 0.06; Fig. 8E), though these observations did not reach our pre-defined significance requirement of P < 0.05 on both tests. Overall, these results indicated that participant 5’s proviral pool was largely stable in terms of genetic composition, though the recently emerging, younger HIV sequences in plasma may have re-seeded the reservoir to a modest extent. No recombinant sequences were recovered from this participant.

Participant 6 acquired HIV around August 2006 and initiated ART in July 2010 (Table 1; Fig. 9A). We inferred 3,000 phylogenies relating 80 plasma env-gp120 sequences isolated at six time points during untreated infection, and 161 proviral env-gp120 sequences collected at 3, 5, 7, and 8 years post-ART. Of these, 2,528 passed QC and yielded a mean root date of January 2007 (95% HPD, from September 2006 to May 2007), which was only slightly later than the clinically estimated root date (Table S1). Proviruses persisting on ART interspersed throughout the whole tree, with the oldest dating to March 2007 (Fig. 9B). We observed a number of long-lived clones, including one that was recovered at every on-ART time point (Fig. 9B, black arrow). Proviral age distributions remained consistent over time on ART (Kruskal-Wallis P = 0.5; Fig. 9D), with no evidence of changing population structure (AMOVA P = 0.7; CC mean P = 0.7; Fig. 9E). Only one recombinant proviral sequence was isolated from participant 6, whose parents dated to the 2 years prior to ART (Fig. S4B).

Plotting the inferred integration dates of all rebound HIV and on-ART proviral sequences on a scaled timeline ranging from estimated transmission to ART-mediated suppression allowed all participants to be visualized together (Fig. 10). These data clearly show that proviruses sampled on ART span a wide age range, which in some cases comprises the entirety of untreated infection (e.g., participants 1, 3, and 4). The figure also shows that rebound viruses, shown as circles for participants 1 and 5, represent a younger subset of the overall persisting proviral pool (Mann-Whitney P = 0.0001).

Though env-gp120 is commonly used for within-host evolutionary studies (6, 9, 10, 40), evolution in this region is characterized by (sometimes dramatic) genetic bottlenecks largely imposed by the evolving antibody response (41). To confirm that our choice of HIV region did not unduly influence our results, we additionally dated participant 1, 3, and 7’s data sets using gag, as pre-ART gag sequences had previously been collected for these individuals (32).

We begin with participant 7, who was not included in the primary analysis because proviruses were only sampled at a single on-ART time point, but whose sequences are nevertheless critical to this validation (Table 1; Fig. S5). We inferred phylogenies relating newly collected env-gp120 and gag on-ART proviral sequences (amplified independently) to the existing pre-ART plasma HIV RNA sequences (32) (see Fig. S5B and C, for representative trees and divergence plots, and Table S1 for summary statistics). The age distributions of distinct proviruses sampled on ART did not significantly differ based on which gene was used for dating (Mann-Whitney P = 0.07; Fig. S5D). Env-gp120 and gag proviral clonal profiles were also highly comparable, even though these regions were amplified separately (Fig. S5E). Of the five gp120 recombinant proviral sequences recovered during ART in participant 7, none contained parents that dated to substantively different eras of infection (Fig. S6). No recombinant gag sequences were identified.

Gag analysis also corroborated the original observations for participants 1 and 3. For participant 1, the percentage clonality based on gag was 23% (versus 33% for env-gp120) while for participant 3, it was 39% (versus 58% for env-gp120). Gag clones also increased over time and were observed across multiple time points (Fig. S7). For participant 1, gag analysis also confirmed that plasma rebound HIV sequences were significantly younger than, and represented a distinctive population compared to, the proviral pool that persisted during ART (Fig. S8E; both AMOVA and CC mean P = 0). Gag analysis also corroborated the independent reactivation of a distinct and slightly more ancestral proviral lineage at the second rebound sampling time point in April 2011 (Fig. S8B; bracket). It also confirmed that persisting proviruses had integrated throughout the entire course of untreated infection and that proviral age distributions were stable over time on ART (Kruskal-Wallis P = 0.6) with no evidence of temporal population structure (AMOVA P = 0.6; CC mean P = 0.62) (Fig. S8D and E). For participant 3, gag analysis similarly confirmed that proviral ages were stable during ART (Kruskal-Wallis P = 0.8) with no evidence of temporal population structure (AMOVA P = 0.74; CC mean P = 0.97) (Fig. S9E).

Importantly, the integration date distributions of on-ART sequences that were inferred from env-gp120 versus gag were highly comparable. The sole exception was participant 1’s initial rebounding population, where env-gp120 analysis returned on average older integration dates than gag analysis (P = 0.005, Fig. S10A). Other than this, env-gp120 and gag analyses produced similar proviral integration dates for all seven proviral time points analyzed for participants 1 and 3 and the remaining two rebound time points analyzed for participant 1 (all P ≥ 0.1, Fig. S10A to C). No recombinant gag sequences were identified in participants 1 and 3. Overall, these observations indicate that our findings are not majorly influenced by the HIV region studied.

As there is uncertainty in phylogenetic reconstruction, we also corroborated our findings using a “tree-free” approach. To do this, we leveraged the shifts in HIV coreceptor usage that occurred over time in participants 1, 2, and 4 (coreceptor shifts did not occur in the other participants). We inferred coreceptor usage from the V3 region of env-gp120 sequences using the geno2pheno (coreceptor) algorithm, which assigns each sequence a false positive rate (FPR) between 0 and 100. The FPR represents the likelihood that a CCR5-using virus is misclassified as CXCR4-using (i.e., sequences with low FPR are more likely to be CXCR4-using). For participants 1 and 2, coreceptor usage shifted from CXCR4-using to CCR5-using during untreated infection (participant 1’s shift was previously documented (32)), while participant 4 had a minority CXCR4-using population in early infection that steadily became more dominant (Fig. S11).

The coreceptor usage results corroborated the phylogenetic ones. Consistent with the tree-based analysis, participant 1’s FPR distributions differed significantly between rebound HIV and proviruses persisting longitudinally on ART (Mann-Whitney P = 0.002 for comparison of all rebound viruses to all proviruses; P = 0.003 for comparison of initial rebound viruses to all proviruses). In contrast, yet also consistent with the tree-based analyses, proviral FPR distributions were stable over time in both participants 1 and 2 (Kruskal-Wallis P = 0.3 (Fig. S11A) and P = 0.2 (Fig. S11B), respectively)) but shifted over time in participant 4 (Kruskal-Wallis P = 0.02, where the average FPR was significantly higher at the latest time point compared with one of the earlier ones; post-test P = 0.04 Fig. S11C).

As one of our objectives was to investigate proviral genetic stability (in terms of distinct lineages) during the initial years of ART, we concluded by testing for trends in overall within-host proviral diversity over time. We used two metrics: the grand mean patristic distance (Fig. 11A) and the mean phylogenetic diversity (Fig. 11B). Both metrics were computed from all distinct, non-recombinant proviral sequences isolated at each time point. For grand mean patristic distance, a Friedman test comparing the first three time points across all six participants yielded P = 0.6, while comparison across all four time points for participants 1, 4, and 6 yielded P = 0.7 (Fig. 11A). The corresponding values for mean phylogenetic diversity were P = 0.4 and P = 0.7 (Fig. 11B). These observations further support the notion that proviral genetic diversity is stable during the initial years of ART.

Of note, the diversity of distinct proviral lineages persisting on ART correlated strongly with the overall plasma HIV diversity generated during untreated infection: when measured as grand mean patristic distance, Spearman’s correlation yielded ρ = 0.86, P = 0.02 (Fig. 11C), and when measured as mean phylogenetic diversity, Spearman’s correlation yielded ρ = 0.82, P = 0.03 (Fig. 11D).

The WIHS was a prospective, multi-center, cohort study of US women living with or without HIV (59–61). WIHS participants were recruited at 10 sites over four time periods starting in 1994. While the WIHS study ended in 2019, most participants continue to be followed under similar protocols as part of the MACS/WIHS Combined Cohort Study (MWCCS) (62). Eligibility criteria and study protocols for the WIHS have been previously described (59–61). Briefly, data were collected using structured in-person interviews and standardized physical and laboratory assessments, with study visits occurring every 6 months. Eligible women had documentation of reactive anti-HIV serology (and if positive a confirmatory test) or, if they were HIV-seronegative, had risk factor(s) for HIV exposure. Baseline sociodemographic characteristics and HIV risk factors were similar between HIV-seropositive and HIV-seronegative women.

We studied seven WIHS participants with the following criteria: documented HIV seroconversion during follow-up, initiated ART during chronic infection, at least four longitudinal pre-ART plasma samples available, and on-ART PBMC samples available. These seven participants represented all WIHS participants who met these criteria. HIV infection dates were calculated as the midpoint between last negative and first seropositive study visits. At the time of the last proviral sampling, participants were a median of 53 (range 49–54) years of age and had been receiving ART for a median of 9 years (range 2.8–12.3 years). In total, we collected longitudinal HIV RNA sequences over a median of 9 (range 2–13) study visits per participant spanning a median of 9 (range 0.8–11.8) years of untreated infection. Plasma HIV RNA sequences for participants 1, 3, and 7 were published previously (32). In addition, we collected single-genome-amplified proviral sequences over a median of 3 (range 1–4) time points per participant spanning a median of 8.7 (range 2.8–12.3) years on ART (only 10 million PBMCs per on-ART time point were available for analysis). We also isolated plasma HIV RNA sequences from participant 1 at three time points after ART was interrupted and from participant 5 when viral control was initially being lost.

Total nucleic acids were extracted from plasma using the NucliSENS EasyMag (BioMerieux, Marcy-l'Étoile, France). If the plasma viral load (pVL) was <2,000 copies/mL, extracts were DNase I-digested (New England Biolabs) to minimize the risk of amplifying proviral DNA. Genomic DNA was extracted from 10 million PBMCs per time point using the QIAamp DNA Mini Kit (Qiagen). Single-genome amplification of a subgenomic HIV region (env-gp120 and gag where applicable) was performed as follows. For plasma HIV RNA, cDNA (generated using NxtScript reverse transcriptase; Roche) was generated using HIV-specific primers and endpoint diluted such that subsequent nested PCR reactions (generated using the Expand High Fidelity PCR system; Roche) yielded no more than 30% positive amplicons (12, 46). Proviral DNA extracts were similarly endpoint diluted and amplified by nested PCR. For env-gp120, first round primers were as follows: 5′-TTAGGCATCTCCTATGGCAGGAAGAAGCGG-3′ (forward; HIV reference strain HXB2 genomic nucleotide coordinates 5957-5986) and 5′-TAAGTCATTGGTCTTAAAGGTACC-3′ (reverse; HXB2 9038-9015); second round primers were 5'- GGCCGCGTCGACAAGAGCAGAAGACAGTGGCAATGA-3′ (forward; HXB2 6194-6228) and 5′-GGCCGCGGATCCGTGCTTCCTGCTGCTCCCAAGAAC-3′ (reverse; HXB2 7823-7787) (32). For gag, first round primers were 5′-AAATCTCTAGCAGTGGCGCCCGAACAG-3′ (forward; HXB2 629-649) and 5′-TAACCCTGCGGGATGTGGTATTCC-3′ (reverse; HXB2 2849-2826); second round primers were 5′-GCAGGACTCGGCTTGCTGAA-3′ (forward; HXB2 691-710) and 5′-TATCATCTGCTCCTGTATC-3′ (reverse; HXB2 2343-2325). Negative controls were included in every run. For extracts from plasma samples with low pVL, we confirmed that HIV RNA amplification did not occur in the absence of reverse transcription. Amplicons were sequenced on a 3730xl automated DNA sequencer using BigDye (v3.1) chemistry (Applied Biosystems). Chromatograms were analyzed using Sequencher (v5.0/v5.4.6) (GeneCodes). Sequences with nucleotide mixtures were excluded from analysis.

Hypermutated sequences were identified using Hypermut (63). Sequences exhibiting evidence of putative within-host recombination were identified using RDP4 v4.1 (38). This program identifies the putative major and minor parent sequences for each recombinant, along with their approximate breakpoints, allowing us to date each sequence fragment by assigning the date of origin of its parent. Hypermutant and recombinant sequences, along with those with minor defects (e.g., small deletions), were retained in the analysis of clonality but excluded from phylogenetic inference. Sequence alignments were performed in a codon-aware manner using MAFFT v7.471 (64). Alignments were inspected and manually edited in AliView v1.26 (65). Following automated model selection using ModelFinder (66), between-host phylogenies were inferred by maximum-likelihood methods using IQ-TREE 2, with the ultrafast bootstrap option (1,000 bootstraps) (67, 68).

Within-host phylogenies were inferred from env-gp120 and gag sequence alignments comprising all plasma and proviral sequences collected per participant. To mitigate uncertainty in phylogenetic reconstruction, we inferred a median of 4,500 (range 1,500–15,000) phylogenies per participant using Bayesian approaches and conditioned results across all trees. To do this, we first reduced each within-host nucleotide sequence alignment to distinct, intact, non-recombinant sequences and determined the best-fitting nucleotide substitution model for each alignment using jModelTest v2.1.10 (Table S1) (69). Next, Markov chain Monte Carlo (MCMC) methods were used to build a distribution of phylogenies per participant without enforcing a molecular clock. Two parallel runs with MCMC chains of a median of 30 million generations, sampled every 10,000 generations, were performed in MrBayes, v3.2.5 (70) using the best-fitting substitution model and model-specific or default priors. Convergence was assessed by ensuring the standard deviation of split frequencies was <0.04, and the effective sampling size of all parameters was ≥200 and by visual inspection of parameter traces using Tracer v1.7.2 (71). In the single case where convergence was not achieved (participant 1, env-gp120), the run was terminated at 100 million generations. The first 25% of MCMC iterations was discarded as burn-in, yielding a minimum of 1,500 (maximum 15,000) env-gp120 and gag phylogenies per participant (Table S1).

We then inferred the integration dates of on-ART sequences of interest using a phylogenetic approach (7). First, each tree was rooted at the location that maximized the correlation between the root-to-tip distances of the pre-ART plasma HIV RNA sequences and their sampling dates, which represents the most recent common ancestor of the data set. We then fit a linear model relating the root-to-tip distances of pre-ART plasma HIV sequences to their collection dates, where the slope of this line represents the average within-host gene-specific evolutionary rate, and the x-intercept represents the root date. Model fit was assessed by comparing the model’s Akaike Information Criterion (AIC) to that of a null model with a zero slope. In order to pass QC, a phylogeny required a ΔAIC ≥ 10 and an inferred root date that was prior to the first plasma sampling. A median of 3,842 (range 1,278–7,218) phylogenies per participant passed QC (Table S1). The linear models from these QC-passed phylogenies were used to convert the root-to-tip distances of on-ART sequences of interest to their integration dates. The custom R script for this method is available at https://github.com/cfe-lab/phylodating. The script was implemented via GNU parallel (72) to run more than one tree at a time. The integration dates were then averaged across all QC-passed phylogenies per participant to produce mean integration date estimates and 95% HPD estimates, computed using R package HDInterval (version 0.2.2).

As phylogenies were inferred using distinct within-host HIV sequences only, identical sequences were then grafted back onto the phylogenies using the add.tips function in the R package phangorn, v2.8.1, as appropriate for each analysis (73). The example phylogeny shown for each participant was the highest likelihood tree among those that passed QC. Phylogenies and highlighter plots were plotted using the R (v4.1.2) package ggtree version 3.21. Node support values were derived from Bayesian posterior probabilities generated from the consensus trees.

Within-host proviral populations sampled during ART were tested for evidence of population genetic structure using AMOVA (74), which is a genetic distance-based test, and the CC test (75), which is a tree-based test. These were chosen because they can test for population structure across more than two time points. Tests were performed on distinct sequences per on-ART time point. AMOVA was implemented in the R package pegas, v1.1 (76) using the K80 substitution model (77), where statistical significance was assessed via 1,000 permutation tests. CC test statistics were averaged over all QC-passed phylogenies per participant, using a custom R script available at https://github.com/brj1/HIVCompartmentalization. A data set was classified as having evidence of a population structure if both tests returned P < 0.05 (where, for CC, this was defined as a mean P < 0.05 over all QC-passed trees).

Within-host proviral diversity was quantified using two metrics: the grand mean patristic (tip-to-tip phylogenetic) distance, calculated as the mean patristic distance between all pairs of distinct sequences per time point averaged over all QC-passed phylogenies, and the mean phylogenetic diversity, calculated by summing the edge lengths of all distinct sequences per time point, averaged over all QC-passed phylogenies (78). Both metrics were computed using custom R scripts available at https://github.com/brj1/HIVCompartmentalization. All tests of significance were two tailed (except compartmentalization tests which are one tailed), with P < 0.05 denoting statistical significance. All other statistical analyses were performed in Prism, v9.0 (GraphPad Software).

Coreceptor usage was predicted from the V3 region of HIV-1 env-gp120 sequences using geno2pheno (coreceptor) (https://coreceptor.geno2pheno.org) (79). This support vector machine-based approach assigns each V3 sequence a FPR that represents the probability of falsely classifying a CCR5-using virus as CXCR4-using. We considered V3 sequences with FPR < 10% as CXCR4-using and those with FPR ≥ 10% as CCR5-using.