Parallel analysis of transcription, integration, and sequence of single HIV-1 proviruses

1Continuously ongoing, high-level viral replication represents the hallmark of untreated HIV-1 infection. Antiretroviral therapy (ART) inhibits active viral replication, but infected cells harboring chromosomally integrated proviruses persist life long and can drive rebound viremia after treatment interruption, necessitating indefinite ART. Such residual HIV-1 proviruses have traditionally been regarded as “transcriptionally silent,” but proviral gene expression is not affected by currently available antiretroviral agents and might be ongoing in some cells. Viral gene expression may represent a key and possibly the only transcriptional feature that distinguishes HIV-1-infected cells from uninfected counterparts; moreover, proviral transcriptional activity can selectively expose HIV-1-infected cells to recognition by adaptive immune responses and increase cell-intrinsic viral cytopathic effects. For these reasons, proviral gene expression might influence the persistence, longevity, and fate of virally infected cells during ongoing suppressive ART.

2Footprints of host immune selection mechanisms can be readily detected in untreated people living with HIV-1 (PLHIV) through the emergence of escape mutations in cytotoxic T cell epitopes and in antibody contact regions. However, after initiation of ART, viral mutational escape is effectively blocked, reducing the probability that selection mechanisms can be inferred from viral sequence variations; instead, selection forces influencing HIV-1-infected cells during suppressive ART may become visible using assays that capture the transcriptional activity and/or the chromosomal location of proviruses. Until recently, HIV-1-infected cells were typically evaluated by PCR assays designed to amplify small segments of HIV-1 DNA, an approach that did not provide information beyond the simple presence or absence of the amplified proviral region. Subsequent technical advances permitted the analysis of individual proviral sequences through single-template, near full-length HIV-1 sequencing; in combination with bioinformatic tools designed to discriminate between genome-intact and defective proviruses, this method allowed for the profiling of proviral landscapes at single-genome resolution, frequently revealing clusters of sequence-identical proviral species derived from clonal proliferation of infected cells. When preceded by whole-genome amplification, single-genome near full-length proviral sequencing can be performed in conjunction with corresponding chromosomal integration site (IS) analysis. This technique, referred to as matched integration site and proviral sequencing (MIP-seq), recently demonstrated a highly disproportionate accumulation of genome-intact proviruses from elite controllers in heterochromatin regions, likely reflecting immune-mediated selection mechanisms that preferentially eliminate proviruses integrated in chromosomal regions more permissive to viral gene expression. These studies indirectly supported the assumption that proviral transcriptional activity may increase the vulnerability of HIV-1-infected cells to host immune activity but failed to provide a direct assessment of the transcriptional behavior of individual proviruses. A simultaneous ex vivo analysis of the transcriptional activity and chromosomal location of proviruses in participant-derived cells would greatly facilitate the understanding of mechanisms that govern proviral gene expression, influence proviral susceptibility to host immune recognition, and determine immune selection of HIV-1-infected cells during ART.

3In this study, we developed an assay for the high-resolution analysis of individual HIV-1-infected cells, designed to simultaneously capture the transcriptional activity, the sequence, and the chromosomal IS of single HIV-1 proviruses. In combination with genome-wide assessments of epigenetic chromatin features from reference data of primary CD4 T cells, this analysis enabled us to measure the frequency of transcriptionally active proviruses, annotate their corresponding chromosomal locations with surrounding epigenetic chromatin features, and profile their longitudinal evolution during continuous ART.

4We adapted a previously described protocol for simultaneous genome and transcriptome analysis (G&T-seq) to develop an assay designed for the parallel interrogation of HIV-1 RNA, the HIV-1 chromosomal IS, and the corresponding proviral sequence of individual HIV-1-infected cells, here termed parallel HIV-1 RNA, integration site, and proviral sequencing (PRIP-seq) (Figures 1A, S1A, and S1B). For this purpose, peripheral blood mononuclear cells (PBMC) purified from persons living with HIV-1 were diluted to single HIV-1-infected cells based on Poisson distribution statistics, lysed, and exposed to biotin-labeled primers designed to anneal to HIV-1 RNA. Subsequently, viral RNA transcripts were physically separated using ultrasensitive immunomagnetic enrichment technologies specifically adjusted for the isolation of low-abundance nucleic acids. Viral RNA was then reverse transcribed, followed by cDNA amplification using a Smart-seq2 protocol. Next, viral cDNA copies were quantified using droplet digital PCR with primers annealing to the HIV-1 long LTR region and to more elongated HIV-1 RNA transcripts, including those containing HIV-1 pol, nef, tat-rev, and poly-A sequences. In parallel, corresponding genomic DNA from each cell lysate was subjected to phi29-catalyzed multiple displacement amplification (MDA), enabling subsequent simultaneous analysis of near full-length proviral sequences and their respective chromosomal IS using a previously described protocol. Briefly, MDA products were split and separately subjected to near full-length HIV-1 amplification as well as amplification of 3′-LTR viral/host junctions by integration site loop amplification (ISLA) assays; resulting amplification products were sequenced on the Illumina MiSeq platform. Technical evaluations of this assay platform with defined copy numbers of HIV-1 RNA transcripts spiked into background populations of HIV-1-uninfected cells revealed that as few as 1–5 viral RNA copies could be detected, and there was a 50% probability of detecting 10 HIV-1 RNA copies. Moreover, there was a strong correlation between input HIV-1 RNA copy numbers and the postamplification cDNA copy numbers, indicating that postamplification viral cDNA copies provide a quantitative estimate of the original number of viral RNA molecules present in HIV-1-infected cells (Figures S1C and S1D). Our protocol also reliably excluded contamination of viral RNA with genomic HIV-1 DNA (Figure S1E).

5We applied this complex experimental workflow to PBMC samples that were longitudinally collected from 5 ART-treated PLHIV (P1–P3, P5, P6); in one study participant (P4), a single cross-sectional PBMC sample collected after long-term ART was analyzed (Figure S2). Collectively, we analyzed a total of 1,270 single HIV-1 proviruses from these study participants, of which 147 encoded genome-intact HIV-1; the remaining proviruses displayed various structural sequence defects precluding active viral replication (Figure 1B; Table S1). The chromosomal IS was successfully identified in a total of 701 of these proviral sequences; failure to determine the proviral chromosomal location in the remaining cells was frequently due to the absence of ISLA primer annealing sites in proviruses with large deletions. IS coordinates were most frequently located in introns of highly expressed genes, consistent with previous results; moreover, among proviruses integrated in genes, opposite orientation relative to the host gene was approximately twice as common as concordant/same orientation (Figures 1C and S1F).

6Among all 1,270 analyzed proviruses, we noted that HIV-1 RNA transcription varied profoundly: HIV-1 long LTR transcripts were detected in 429 (33.8%) of all 1,270 analyzed HIV-1 proviruses when all sequences were considered, and in 296 (31.1%) of 953 analyzed proviruses when sequences from aviremic time points were selectively evaluated (Figure 1D). Because structural sequence defects in the viral promoter region may compromise the ability of proviruses to express HIV-1 RNA, we used MDA products to amplify and sequence the viral core promoter region in a subset of HIV-1-infected cells; restricting the analysis to proviruses with confirmed intact promoter regions did not markedly affect the relative proportion of transcriptionally active proviruses (Figure 1D). To explore the transcriptional activity of proviruses in the context of epigenetic and architectural features of human genome organization, we annotated all IS coordinates with genome-wide data derived from assays to characterize host gene expression (RNA-seq), chromatin accessibility (ATAC-seq), and three-dimensional (3D) chromosomal contacts (in situ Hi-C). RNA-seq and ATAC-seq data were previously generated from primary total CD4 T cells of three ART-treated PLHIV; corresponding Hi-C data from the same three study participants were generated using an in situ Hi-C protocol described previously (Figure S1G). In addition, we used ChIP-seq data corresponding to activating and inhibitory histone modifications from primary total CD4 T cells included in the ROADMAP database; genome-wide cytosine methylation data from the iMethyl database were also evaluated to determine epigenetic cytosine methylation in defined regions upstream of HIV-1 promoters. Collectively, these data permitted a comprehensive investigation of the sequence and transcriptional activity of individual proviruses, combined with an analysis of their corresponding chromosomal locations relative to a diverse range of transcriptional, genomic, and epigenetic host features.

7We initially investigated the transcriptional behavior of HIV-1 proviruses located in nongenic regions, which are typically disfavored for proviral integration because of host factors that bias chromosomal insertion sites toward active transcription units in the human genome. Among all 701 proviruses for which chromosomal IS were identified, 117 (17%) were located in such nongenic/pseudogenic locations (Figure 1E); intact proviruses were more frequently located in nongenic regions relative to defective proviruses (32/122 intact versus 85/579 defective, p = 0.003). Notably, HIV-1 proviruses integrated in intergenic regions of the human genome displayed significantly lower transcriptional activity (Figures 1E–1G); this was particularly true for the small number (n = 34) of HIV-1 proviruses integrated in repetitive satellite or microsatellite DNA, of which the vast majority (94%) were completely transcriptionally silent, consistent with a state of “deep viral latency”. Because intergenic DNA is not generally resistant to transcription and RNA complementary to nongenic and satellite DNA is detectable in many species, we explored mechanisms that might explain the repressed transcriptional activity of proviruses integrated in nongenic DNA. As expected, we noted that RNA-seq reads surrounding nongenic locations of integrated HIV-1 DNA were significantly diminished, relative to genic sites harboring integrated HIV-1 DNA (Figure S3A). A similar finding was made for ChIP-seq reads corresponding to activating histone features (H3K27ac, H3K4me1, and H3K4me3); however, no difference was noted for ChIP-seq reads related to inhibitory histone marks (H3K27me3 and H3K9me3) (Figures S3B and S3C). Compared with those in genes, nongenic proviral IS were also in regions with significantly reduced chromatin accessibility, as determined by ATAC-seq (Figure S3D). Interestingly, an alignment of HIV-1 IS to spatial 3D chromatin organization features, evaluated by genome-wide in situ Hi-C, demonstrated that nongenic chromosomal regions harboring HIV-1 IS displayed significantly increased chromosomal distances to frequently interacting regions (FIREs) and to topologically associated domains (TADs) (Figures S3E and S3F); moreover, the nongenic regions containing integrated HIV-1 DNA had reduced numbers of 3D intra- and interchromosomal chromatin contacts (Figures S3G–S3H). Together, these results demonstrate that proviruses in nongenic locations have significantly weaker viral transcriptional activity, likely because of nonpermissive genomic and epigenetic chromatin features at chromosomal IS.

8To dissect global features associated with the transcriptional activity of proviruses, we analyzed epigenetic characteristics in DNA regions linearly surrounding IS in genic and nongenic regions (Figures 2A–2D and S4A–S4D). Although positioning of IS relative to the most proximal host transcriptional start sites (TSSs) had no noticeable influence on proviral transcription (Figures S4A and S4B), we observed that ATAC-seq reads and ChIP-seq reads corresponding to activating histone modifications (H3K4me1, H3K4me3, and H3K27ac) in linear proximity of HIV-1 DNA IS were significantly higher for transcriptionally active proviruses (Figures 2B, 2C, and S4C); this was true when all proviruses from all time points were considered and when proviruses in genic locations or from aviremic time points were selectively analyzed. Proviral transcriptional silence, on the other hand, was associated with hypermethylated cytosine residues in genomic DNA upstream of the 5′-LTR junction of proviral IS (Figure 2D); inhibitory H3K27me3-specific ChIP-seq marks appeared unrelated to proviral transcriptional behavior (Figure S4D). In a subsequent analysis, we investigated proviral transcriptional activity relative to genomic and epigenetic features in intra- and interchromosomal 3D contact regions, determined by in situ Hi-C (Figures 2E–2I). The total number of significant 3D contacts tended to be higher in regions harboring IS of transcriptionally active proviruses, and their chromosomal distances to FIREs were reduced, suggesting that proviral gene expression is facilitated by integration into more interactive chromosomal regions (Figures 2F, 2G, S4E, and S4F). Moreover, there were trends for higher levels of activating histone-specific ChIP-seq reads (Figures 2H and S4H) and ATAC/RNA-seq reads (Figures 2I and S4G) in 3D contact regions of IS of transcriptionally active proviruses. Hypothesizing that combinations of chromatin features in linear, intrachromosomal, and interchromosomal contact regions of chromosomal IS may influence proviral transcriptional behavior, we calculated the sum of reads corresponding to epigenetic features in all three of these chromatin compartments for each individual proviral IS; this composite marker may better approximate the totality of epigenetic signals influencing a given provirus in the 3D spatial configuration of the genome. Collectively, these analyses supported the assumption that combined epigenetic features in linear and 3D contact regions of IS can jointly regulate the transcriptional activity of HIV-1 DNA (Figures 2J, 2K, and S4I). We also observed some evidence for privileged interchromosomal interactions between genomic positions of transcriptionally active proviruses, possibly reflecting the localization of transcriptionally active proviral species at highly permissive and interactive chromatin regions on the outer surface of chromosomal territories (“transcription factories”) (Figure S4J); however, signs of such interchromosomal “transcriptional interactomes” were not observed genome wide, possibly because of the disproportionate overrepresentation of the small, gene-rich chromosomes in 3D interchromosomal interactions. Collectively, these data suggest that the relative position of proviruses in the nuclear 3D chromatin architecture can influence viral transcription through combined ensemble effects of cis and trans epigenetic signals and that HIV-1 can take advantage of the full complexity of the human genome and epigenome in its 3D spatial configuration to fine-tune viral transcriptional behavior.

9Together, these results identify several genomic and epigenetic features permitting discrimination of transcriptionally active proviruses across a large pool of proviral species isolated from a heterogeneous group of participants and sampling time points; however, the discriminatory effect size of each of these parameters, when assessed as isolated variables, was relatively weak. We subsequently considered a combinatorial analysis of genomic and epigenetic parameters for which significant associations with transcriptional activity were noted in univariable comparisons; these investigations indicated a progressive increase in transcriptional activity among proviruses stratified according to combinations of distinct IS features (Figure 2L). Consistent with this observation, a formal multivariable logistic regression model demonstrated that a combinatorial analysis of genomic and epigenetic IS characteristics facilitated discrimination of transcriptionally active proviruses (Figures 2M–2O) and supported positive associations between activating epigenetic chromatin features and proviral transcriptional activity. However, effects in this model varied considerably with the study subject and were critically influenced by the sampling time point, with longer durations of ART being generally associated with more limited proviral gene expression (Figure 2O). Therefore, the regulation of proviral transcriptional behavior may follow distinct evolutionary pathways in individual subjects, requiring participant-specific longitudinal analyses for further clarification.

10To investigate the time-dependent evolution of proviruses during suppressive ART, we evaluated the trajectory of HIV-1-infected cells around the time of ART initiation and during the ensuing 9–12 years of therapy in P1 and P2; in P3, data from year 1 and year 9 after ART commencement were available for longitudinal investigation. The frequency of all proviruses tended to decline after ART initiation in these three participants (Figures S5A–S5L). However, this decrease was substantially more pronounced for transcriptionally active proviruses, categorized by the detection of any HIV-1 RNA transcript, high-level HIV-1 RNA expression (>10,000 postamplification copies), or elongated HIV-1 RNA transcripts (Figures 3A–3C, S5A, S5E, and S5I). The accumulation of transcriptionally silent proviruses over time coincided with a gradual increase in the frequency of proviral IS in nongenic and satellite DNA (Figure 3D); however, the progressive proportional reduction in HIV-1 RNA-expressing cells was similarly visible among proviruses integrated in genic chromosomal regions (Figures 3B, S5B, S5F, and S5J) and resulted in profound alterations in the composition of the proviral landscape, with a disproportionate overrepresentation of transcriptionally silent proviruses after prolonged ART (Figure 3E). This pattern was consistent with low proportions of transcriptionally active proviruses detected cross-sectionally in study person 4 after 15 years of ART (Figures S5M–S5O). Notably, the longitudinal decline of HIV-1 RNA-expressing proviruses was more obvious for intact HIV-1 DNA compared with defective proviral species, leading to a markedly biased proviral IS profile with almost completely undetectable transcriptionally active intact proviruses after long-term ART (Figures 3F–3H). By contrast, the contribution of transcriptionally silent intact as well as transcriptionally silent defective proviruses to the total pool of HIV-1-infected cells remained relatively stable or expanded over the entire observation period (Figure 3H). Large clones of proviral sequences were detected in two of these study participants after several years of continuous ART. In P2, we noted a large genome-intact clone integrated in the ZNF140 gene and a second clone of near full-length proviruses with a premature stop codon (PSC) in gag integrated in pericentromeric satellite DNA on chromosome 16; both of these proviral clones showed no or minimal transcriptional activity (Figure 5). In P3, we noticed the evolution of two clones of intact proviruses, integrated in centromeric satellite DNA on chromosome 10 and on chromosome Y, respectively; no members of these two clones expressed detectable HIV-1 RNA (Figure 5). Together, these observations suggest a progressive selection of proviruses with lower transcriptional activity and chromosomal integration into nongenic or satellite DNA locations during prolonged ART.

11In contrast to P1–P3, study participant 5 displayed a biphasic evolution of HIV-1-infected cells, with an initial decline occurring during the first year of ART and a subsequent increase during prolonged ART (Figures 3I–3K). A closer investigation demonstrated that this longitudinal pattern was primarily driven by large clones of virally infected cells (Figures 3K–3M) that became visible after prolonged ART durations and frequently expressed viral long LTR HIV-1 RNA and, to a lesser extent, elongated HIV-1 RNA transcripts; in at least two cases, these large, transcriptionally active clones of cells encoded genome-intact HIV-1 (Figure 5). IS of these large, transcriptionally active proviruses were located at chromosomal positions with strong activating epigenetic support (Figure 6). By contrast, among proviruses with IS detected only once (derived from smaller clones or nonclonal HIV-1-infected cells), the proportion of proviruses expressing HIV-1 RNA declined over time (Figure 3J); this longitudinal decrease was most visible in the context of elongated HIV-1 RNA transcripts. A similar pattern was noted in P6, who had detectable low-level plasma viremia (>20 but <100 copies/mL) despite prolonged treatment with ART, in the absence of notable proviral resistance mutations against antiretroviral drugs. In this person, we generally noted a high (>50%) proportion of transcriptionally active proviruses (Figures 3N–3P). Nevertheless, there was a trend for a longitudinal decline of proviruses that were detected once and expressed elongated and/or high-level viral RNA (Figure 3O). By contrast, no such decrease was observed for large proviral clones that remained transcriptionally active and were integrated in chromosomal locations with strong activating chromatin features (Figures 3P–3R, 5, and 6); two of these transcriptionally active clones encoded genome-intact HIV-1 DNA (Figure 5). Collectively, the longitudinal trajectories in P5 and P6 resulted in a proviral landscape characterized by a gradual decline of HIV-1 RNA-expressing proviruses originating from small clones or nonclonal HIV-1-infected cells, coupled with a parallel expansion of large proviral clones frequently displaying strong proviral transcriptional activity. These results are generally consistent with a progressive selection advantage for proviruses with lower transcriptional activity during extended ART. However, large transcriptionally active proviral clones in P5 and P6 appeared to violate this evolutionary pattern and persisted/expanded despite strong expression of viral RNA, possibly because they can outcompete negative host selection forces through enhanced cell proliferation.

12To better understand longitudinal selection mechanisms underlying changes in proviral transcriptional activity, we analyzed proviral chromosomal IS features over time in P1, P2, and P5, in whom PRIP-seq data were obtained at relatively similar time points (viremia, 1 year after ART initiation, and 9–12 years after ART initiation). Aligning IS coordinates to genome-wide cytosine methylation data in reference CD4 T cells, we noted a progressive longitudinal accumulation of genic proviruses integrated in chromosomal regions with hypermethylated cytosine residues in host DNA upstream of the proviral 5′-LTR promoter, suggesting a role of epigenetic cytosine methylation in silencing proviral transcriptional activity during prolonged ART (Figures 4A and 4B). Although progressive transcriptional silencing was similarly visible for genic proviruses in same and opposite configurations to the most proximal host TSS (Figures S5P and S5Q), we noted marked differences in the longitudinal evolution of IS features for proviruses in the same versus opposite orientations to the nearest host TSS. Although chromosomal distances between proviral IS and the nearest same-oriented TSSs tended to increase over time, this trend was less apparent for chromosomal distances to the nearest opposite-oriented TSSs (Figure 4C). Moreover, using RNA-seq data from primary CD4 T cells from reference datasets, we noted that the host gene expression intensity at the nearest same-oriented TSSs remained stable over time, whereas host transcriptional activity at the nearest opposite-oriented TSSs increased, most notably in the context of TSSs in convergent orientation to the provirus (Figure 4D). In addition, among proviruses with the same orientation as the nearest host TSS, ATAC-seq and activating ChIP-seq reads in linear proximity of regions harboring integrated HIV-1 DNA decreased significantly over time (Figures 4E–4G); by contrast, there was a trend for increasing levels of ATAC-seq and activating ChIP-seq reads in linear proximity of IS for proviruses in opposite orientation to the nearest host TSS (Figures 4H–4J). Together, these data suggest that progressive selection of more transcriptionally silent proviruses over time is partly achieved by two complementary mechanisms: (1) selection forces that promote preferential persistence of proviruses with relatively increased distance to same-oriented host TSSs; transcriptional silencing of these proviruses is likely due to deprived proviral access to activating epigenetic chromatin signals and to host transcriptional machinery; and (2) host factors that select for proviruses in relative proximity to highly expressed host TSSs in opposite orientation; proviral transcriptional repression in this context may be achieved through transcriptional interference from host gene expression.

13Clusters of clonal, sequence-identical HIV-1 proviruses were detected in all six study subjects, although their frequencies varied considerably; these sequence clusters (n = 42 in total) result from clonal proliferation of virally infected cells (Figure 5). Particularly, we noted multiple large clones integrated in nongenic/pseudogenic DNA that evolved after prolonged ART and displayed no or minimal signs of transcriptional activity, specifically in P1–P4. Such clusters of sequence-identical proviruses included a clone with a large deletion integrated in centromeric satellite DNA (P1), a clone with a PSC in pericentromeric satellite DNA (P2), two genome-intact clones in centromeric satellite DNA (P3), and a clone with a 5′-LTR defect in a pseudogenic region (P4). Moreover, one large genome-intact clone integrated in the ZNF140 gene on chromosome 4 was observed after 12 years of treatment in P2; this clone was also completely transcriptionally silent across all 8 member sequences (Figure 5). By contrast, proviral clones with transcriptionally active members in P5 and P6 were frequently integrated in chromosomal locations surrounded by exceptionally strong activating epigenetic chromatin signals in their immediate chromatin proximity, particularly with regard to H3K4me1-specific ChIP-seq reads (Figures 5 and 6A–6E). Similarly, high levels of RNA-seq, ATAC-seq, and activating ChIP-seq signals were observed for transcriptionally active clones when sums of epigenetic sequencing reads in linear proximity and in 3D interchromosomal contact regions of IS were evaluated as composite readouts (Figures 6F–6H). Notably, one large transcriptionally active intact clone was detected in a nongenic region in P6; the IS of this clone was again located in immediate proximity to activating histone modifications, although there were no detectable RNA-seq reads at this nongenic site (Figures 5 and 6B).

14The persistence and expansion of large, transcriptionally active proviral clones during long-term ART appeared to violate the longitudinal trajectory of proviral sequences detected once, for which we documented a progressive accumulation of transcriptionally silent proviruses with features of deeper transcriptional latency over time; this raises the question of how transcriptionally active proviruses can resist negative host selection mechanisms. Notably, genes harboring transcriptionally active clonal proviruses were frequently involved in the regulation of cell proliferation and oncogenesis, consistent with previous findings. This was particularly true for a large, transcriptionally active proviral clone integrated into the KDM2A gene (Figures 5 and 6A), which encodes a lysine-specific demethylase that can enhance cell-autonomous proliferation through downregulation of TET2; disruption of TET2 has previously been associated with enhanced clonal proliferation of CAR-encoding T cells. An oncogenic role was also reported for other genes harboring transcriptionally active proviral clones, such as CFAP299/C4orf22, USP34, and GON4L (Figures 5 and 6A). Retroviral integration into such cancer/proliferation-associated genes may result in splicing-induced chimeric host/viral RNA, leading to fusion proteins that can stimulate cell-autonomous clonal proliferation. For other transcriptionally active clonal proviruses, oncogenic functions of genes harboring the proviral IS were less obvious, and proliferation of the corresponding HIV-1-infected cells may be driven by antigen-specific effects. Independently of the mechanisms that drive proliferation of clonal HIV-1-infected cells, we propose that transcriptionally active proviral clones can persist long-term during ART when negative selection forces are outperformed by elevated cell turnover rates that effectively replenish and expand the clonal pool size.

15Latency-reversing agents (LRAs) have been developed for the purpose of enhancing proviral transcriptional activity, with the intention of increasing their susceptibility to host immune factors; however, the influence of such LRAs on transcriptional activity may vary among individual proviruses and has rarely been assessed in single proviruses. To address this, we used the PRIP-seq platform to evaluate the transcriptional activity of HIV-1-infected cells in the presence or absence of stimulation with phorbol 12-myristate 13-acetate (PMA) and ionomycin, activators of protein kinase C (PKC) that are frequently considered the most powerful inducers of HIV-1 gene expression and commonly used as the positive control in HIV-1 reactivation assays. Briefly, participant-derived CD4 T cells were diluted to single HIV-1-infected cells, stimulated with PMA/ionomycin or with control medium for 12 h, and subsequently subjected to PRIP-seq assays. Overall, we assayed a total of 222 proviruses from two different ART-treated participants (P4 and P7); 97 cells were analyzed after stimulation with PMA and 125 cells after treatment with medium alone (Table S2). We observed no significant differences between the proportion of cells producing detectable levels of HIV-1 long LTR or elongated transcripts in stimulated versus nonstimulated cells (Figure 7A); however, among transcriptionally active proviruses, the per-cell levels of any HIV-1 RNA and of elongated HIV-1 RNA were higher after stimulation (Figures 7B and 7C). Remarkably, though, 79.4% of analyzed HIV-1-infected cells failed to produce detectable HIV-1 RNA despite stimulation with PMA/ionomycin; such stimulation-refractory proviruses had IS located in closer proximity to ChIP-seq peaks related to the inhibitory histone modifications H3K9me3 and H3K27me3 (Figure 7D), suggesting that only a subset of proviruses with favorable epigenetic chromosomal IS features can respond to PMA/ionomycin. Notably, for 4 clusters of clonal proviruses, we detected at least one proviral species in both stimulated and unstimulated conditions, allowing for a direct comparison of the transcriptional activity of proviruses with shared IS in the presence or absence of stimulation with PMA/ionomycin (Figures 7E and 7F). In this paired analysis, we observed that proviruses with detectable HIV-1 RNA at baseline uniformly showed higher per-cell levels of viral transcripts after stimulation (Figure 7F). By contrast, a member of a proviral clone with an intact HIV-1 promoter integrated in a pseudogenic position (Chr5: 177901756) that did not produce HIV-1 RNA in the absence of stimulation also remained transcriptionally silent after stimulation (Figure 7F); 5 transcriptionally silent members of this clone were also detected in direct ex vivo assessments (P4, Figure 5). Together, these results support the hypothesis that chromosomal location critically determines the susceptibility of proviruses to LRAs; PMA/ionomycin stimulation appears to enhance per-cell levels of viral RNA transcripts but is unable to effectively overcome epigenetic transcriptional blocks and “deep” viral latency in proviruses that are transcriptionally silent at baseline. Future studies will be necessary to determine viral reactivation after stimulation with alternative LRAs.

16Although HIV-1 integration in chromosomal DNA is heavily biased toward active host transcription units, the virus can integrate almost anywhere in the human genome, leading to a diverse portfolio of proviral species that differ profoundly with regard to chromosomal location and transcriptional activity. Seeding such a variegated reservoir of proviral sequences that exist in distinct transcriptional activation states may offer the highest flexibility to resist host immune activity and can be viewed as a probabilistic bet-hedging strategy designed to maximize viral persistence in unpredictable environments. However, how proviral transcriptional behavior and chromosomal location synergize in determining the evolutionary fate and the persistence or elimination of HIV-1-infected cells during ART is unclear. Here, we used an assay to simultaneously profile the transcriptional behavior, the proviral sequence, and the chromosomal position of HIV-1 proviruses in unmanipulated participant-derived cells, permitting the global mapping of transcriptionally active and silent proviruses relative to genome-wide genomic and epigenetic chromatin features and the identification of proviral species that are under evolutionary selection pressure during ART in vivo. We note that instead of being “stable” and “transcriptionally silent,” viral reservoir cells are frequently transcriptionally active and dynamically evolving in response to host-dependent selection forces.

17Remarkably, our work suggests positive selection of proviruses with lower transcriptional activity during prolonged ART, arguably due to preferential elimination of transcriptionally active proviruses with higher vulnerability to host immune activity. These footprints of host selection strongly suggest that HIV-1-infected cells are subject to host immune surveillance during ART, likely because proviral gene expression can be sensed by host immune factors. In contrast to immune selection effects during untreated HIV-1 infection that can be readily detected by accumulation of viral sequence variations, the selection footprints reported here are indirectly inferred through complex profiling of individual HIV-1-infected cells; identifying the specific immune responses that drive such selection processes will represent a major challenge for future studies. Of note, we repeatedly noticed that large proviral clones can resist evolutionary selection processes that otherwise promote longitudinal accumulation of proviral sequences with progressively decreasing transcriptional activity. Although long-term persistence of transcriptionally active proviruses has also been suggested by other investigators, our results suggest that they represent an exception of the rule and are facilitated by two coinciding factors: the presence of strong activating epigenetic signals in linear and 3D contact regions of their IS and an elevated rate of cell turnover that can antagonize negative host selection mechanisms through compensatory repopulation and expansion of the clonal pool size. In addition, upregulation of antiapoptosis and cell survival markers during clonal proliferation of T cells may protect large cell clones against antiviral host immune responses. It is possible that these intact, transcriptionally active proviruses play a critical role for driving viral rebound in case of treatment interruptions; future studies will be necessary to address this.

18Manipulating the transcriptional behavior of proviruses to enhance their susceptibility to immune-mediated elimination, colloquially referred to as the “shock and kill” approach, has been evaluated in multiple clinical trials. Our data suggest that pharmacological increases of proviral transcriptional activity may augment proviral vulnerability to host selection mechanisms, but such effects may have remained invisible in previous clinical trials that primarily used quantitative viral reservoir measurements in bulk CD4 T cells as study endpoints. Indeed, it remains possible that LRAs or immunotherapeutic interventions intensify or accelerate naturally occurring immune selection of proviruses and promote a reservoir structure characterized by features of deeper latency and lower transcriptional activity. These more subtle selection effects may result in a proviral reservoir landscape that is more easily controlled by antiviral host immunity and less likely to fuel rebound viremia; ultimately, such a modulated reservoir structure may allow for a peaceful coexistence between HIV-1 and the host and enable drug-free control of HIV-1 infection. Future experiments on the PRIP-seq platform may help to better dissect the effects of LRAs on single HIV-1-infected cells in samples from clinical trial participants and to detect specific susceptibilities and vulnerabilities of HIV-1 reservoir cells.

19A limitation of our work is that transcriptional profiles of individual proviruses were related to chromatin features from reference datasets, and not from single-cell assessments in infected cells. Technical progress in multidimensional single-cell profiling may in the future permit to address this point; in fact, recent studies have started to evaluate lentiviral chromosomal IS coordinates within ATAC-seq reads from single infected cells. In addition, the transcriptional activity of proviruses is known to be burst driven, and it is likely that transcriptional noise may have influenced our findings, although our work identifies clear contributions of the genomic environment to the transcriptional regulation of HIV-1 proviruses.

20Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Mathias Lichterfeld (mlichterfeld@partners.org).

21This study did not generate new unique reagents.

22HIV-1-infected study participants were recruited at the Massachusetts General Hospital in Boston, MA. PBMC samples were obtained according to protocols approved by the Institutional Review Board. Clinical characteristics of study participants are summarized in Figure S2. Demographic characteristics of the study patients are as follows:

23PBMC isolated according to standard protocols were subjected to DNA extraction using commercial kits (Qiagen DNeasy Blood and Tissue Kit, #69504). We amplified total HIV-1 DNA using droplet digital PCR (Bio-Rad), using primers and probes described previously (127 bp 5’-LTR-gag amplicon; HXB2 coordinates 684-810). PCR was performed using the following program: 95°C for 10 min, 45 cycles of 94°C for 30s and 60°C for 1 min, 98°C for 10 min. The droplets were subsequently read by a QX200 droplet reader and data were analyzed using QuantaSoft software (Bio-Rad).

24Cryopreserved PBMCs were thawed, washed, diluted in 96-well plates to single HIV-1-infected cells according to ddPCR results, so that one virally-infected cell was present in approximately 30% of wells, and lysed with Buffer RLT Plus (Qiagen, #1053393). Subsequently, cell lysates were incubated on a thermal mixer (Eppendorf ThermoMixer C, #5382000023) at 1200 rpm for 20 min at room temperature with complexes of magnetic streptavidin beads (Invitrogen Dynabeads MyOne Streptavidin C1, #65002) linked to biotinylated primers targeting defined regions (poly-A, tat-rev, nef, pol, long LTR) of HIV-1 RNA; primer sequences were previously described and modified to permit subsequent amplification with a Smart-seq2 protocol. Primer sequences are listed in Table S3. After primer annealing, viral RNA was magnetically separated from genomic DNA according to the G&T-Seq protocol described previously. Briefly, the magnetic beads bound to HIV-1 RNA were immobilized to the bottom of each well with an external magnet and the supernatant containing the genomic DNA fraction was transferred to a new plate. The beads were washed twice in a buffer containing 50 mM Tris-HCl, 7.5 mM KCl, 3 mM MgCl2, 10 mM DTT and 0.5% Tween 20; after each wash, the supernatants from each well were collected and added to the corresponding genomic DNA solutions. In a modification of the RNA-binding process for participants 2, 3, 5 and 6, biotinylated poly-A, nef, pol and long LTR primers were added sequentially in 10-minute intervals, rather than simultaneously, since this facilitated the capture and later detection of elongated HIV-1 RNA species. In selected cases, these experiments were performed using an RNA standard (ATCC, VR-3245SD) with known viral copy numbers spiked into a background population of 10,000 cells from an HIV-1-negative person for technical validation of the assay workflow.

25The bead-bound HIV-1 RNA was subjected to reverse transcription with SuperScript II Reverse Transcriptase (Invitrogen, #18064022) using a master mix containing 10 mM dNTPs, 6 mM MgCl2, 1 M betaine, 5 mM DTT, 1 U/μL SUPERase⋅In RNase Inhibitor (Invitrogen, #AM2694), 2 μM template-switching oligonucleotide (5’-AAGCAGTGGTATCAACGCAGAGTACATrGrG+G-3’) and 20 U/μL reverse transcriptase in SuperScript II RT First-Strand Buffer. This mix was incubated on a ThermoMixer at 1500 rpm for 60 min at 42 °C, 30 min at 50 °C and 10 min at 60 °C. Subsequently, cDNA was amplified by PCR with KAPA HiFi HotStart ReadyMix (Roche, #7958935001) and 300 nM ISPCR primer (5’-AAGCAGTGGTATCAACGCAGAGT-3’) using a modified Smart-seq2 protocol, as described in the G&T-Seq procedure; for participants 2, 3, 5 and 6, HIV-1-specific primers for poly-A, nef, tat-rev, pol and long LTR viral transcripts were added after the initial 10 Smart-seq2 amplification cycles to complete an additional 25 amplification cycles. Amplified cDNA was subjected to absolute quantification using ddPCR with primers and probes targeting different regions of HIV-1 transcripts, as previously described. All primer sequences are listed in Table S3.

26Isolated gDNA was incubated with AMPure XP magnetic beads (Beckman Coulter, #A63882) and immobilized at the side of each well with an external magnet. After several washes with 80% ethanol, each well was subjected to multiple displacement amplification (MDA) with phi29 polymerase (Qiagen REPLI-g Single Cell Kit, #150345 or REPLI-g Advanced Single Cell Kit, #150365) for 4 hours (Kit #150345) or for 2 hours (Kit #150365), per the manufacturer’s protocol for amplification of genomic DNA from single eukaryotic cells. After this unbiased whole genome amplification, DNA was again incubated with AMPure XP beads and washed and eluted in water. gDNA from each well was split and separately subjected to viral sequencing and integration site analysis, as described below.

27DNA resulting from whole-genome amplification reactions was subjected to HIV-1 near full-genome amplification using a near full-length amplicon as well as a nonmultiplexed overlapping 5-amplicon approach, as described before. Additionally, amplification of the HIV-1 core promoter and enhancer region was performed, using one of two nested PCRs: one with primer sequences corresponding to HXB2 coordinates 350-372 and 642-661 (first round) and 367-385 and 626-643 (second round); or another nested PCR with primer sequences corresponding to HXB2 coordinates 24-50 and 936-962 (first round) and 76-100 and 797-818 (second round). All primer sequences are listed in Table S3. These promoter PCRs were performed with the following program: 94 °C for 2 min, 30 cycles of 94 °C for 15 s, 60 °C for 30 s, 68 °C for 1 min. PCR products were visualized by agarose gel electrophoresis (Quantify One and ChemiDoc MP Image Lab, BioRad). All near full-length or 5-amplicon positive PCR products, as well as a large number of proviruses with evident major deletions, were subjected to Illumina MiSeq sequencing at the MGH DNA Core facility. Resulting short reads were de novo assembled using Ultracycler v1.0 and aligned to HXB2 to identify large deleterious deletions (<8000bp of the amplicon aligned to HXB2), out-of-frame indels, premature/lethal stop codons, internal inversions, or packaging signal deletions (≥15 bp insertions and/or deletions relative to HXB2), using an automated in-house pipeline written in Python programming language (https://github.com/BWH-Lichterfeld-Lab/Intactness-Pipeline), consistent with prior studies. Presence/absence of APOBEC-3G/3F-associated hypermutations was determined using the Los Alamos National Laboratory (LANL) HIV Sequence Database Hypermut 2.0 program; hypermutated proviruses that also contained other structural defects were categorized as hypermutated when binning proviruses according to proviral genomic defects. Viral sequences that lacked all mutations listed above were classified as “genome-intact” sequences. Proviral sequences identified as “genome-intact” with this algorithm were also classified as “genome-intact” using an alternative analysis procedure. Sequence alignments were performed using MUSCLE and Geneious Prime 2021.0.3 (geneious.com). Phylogenetic distances between sequences were examined using UPGMA trees in Geneious Prime. Viral sequences were considered clonal if they had identical integration sites and 3 or fewer mismatches between proviral sequences; single nucleotide variations in primer binding sites were not considered for clonality analysis.

28Integration sites associated with each viral sequence were obtained using integration site loop amplification (ISLA), using a protocol previously described and DNA produced by whole-genome amplification as template. Resulting PCR products were subjected to next-generation sequencing using Illumina MiSeq. MiSeq paired-end FASTQ files were demultiplexed; small reads (142 bp) were then aligned simultaneously to human reference genome GRCh38 and HIV-1 reference genome HXB2 using bwa-mem. Biocomputational identification of integration sites was performed according to previously-described procedures: Briefly, chimeric reads containing both human and HIV-1 sequences were evaluated for mapping quality based on (i) HIV-1 coordinates mapping to the terminal nucleotides of the viral genome, (ii) absolute counts of chimeric reads, and (iii) depth of sequencing coverage in the host genome adjacent to the viral integration site. The final list of integration sites and corresponding chromosomal annotations was obtained using Ensembl (v86, www.ensembl.org), the UCSC Genome Browser (www.genome.ucsc.edu) and GENCODE (v25, www.gencodegenes.org). Repetitive genomic sequences harboring HIV-1 integration sites were identified using RepeatMasker (www.repeatmasker.org). Cells in which multiple HIV-1 integration sites were detected, either as a result of multiple HIV-1-infected cells being present in one limiting-dilution well or multiple HIV-1 proviruses being integrated into a single cell, were excluded from analysis; the latter possibility is unlikely.

29To evaluate HIV-1 RNA expression in response to stimulation, we isolated CD4+ T cells from participant-derived PBMC (Stemcell EasySep™ Human CD4+ T Cell Isolation Kit, #17952), followed by dilution to single HIV-1-infected cells for incubation with 81 nM PMA and 1.34 μM ionomycin (Biolegend Cell Activation Cocktail (without Brefeldin A), #423302) with 20 U/μL recombinant IL-2 and 400 nM AZT; simultaneously, unstimulated autologous cells were plated in the same way for control purposes. After 12 hours of incubation, stimulated and unstimulated CD4+ T cells were subjected to the PRIP-Seq protocol for analysis of HIV-1 integration sites, proviral sequences and HIV-1 RNA expression from individual HIV-1-infected cells.

30PBMC from three ART-treated HIV-1 participants were used for parallel analysis of CD4 T cells by RNA-Seq, ATAC-Seq, and Hi-C, as described below. ChIP-Seq data were obtained from primary memory CD4 T cells included in the ROADMAP database or from activated CD4 T cells included in the ENCODE project. Cytosine methylation data were obtained from the iMethyl database.

31Briefly, total PBMC were stained with monoclonal antibodies to CD4 (clone RPA-T4, Biolegend, #300518), CD3 (clone OKT3, Biolegend, #317332), CD45RO (clone UCHL1, Biolegend, #304236) and CCR7 (clone G043H7, Biolegend, #353216). Afterwards, cells were washed and CD45RO+ CCR7+ (central-memory) and CD45RO+ CCR7- (effector-memory) and CD3+ CD4+ (total) CD4+ T cells were sorted in a specifically designated biosafety cabinet (Baker Hood), using a FACS Aria cell sorter (BD Biosciences) at 70 pounds per square inch. Cell sorting was performed by the Ragon Institute Imaging Core Facility at MGH and resulted in isolation of lymphocytes with the defined phenotypic characteristics of >95% purity.

32RNA was extracted from total CD4+ T-cell populations using a PicoPure RNA Isolation Kit (Applied Biosystems, #0204). RNA-Seq libraries were generated as previously described. Briefly, whole transcriptome amplification and tagmentation-based library preparation was performed using Smart-seq2, followed by sequencing on a NextSeq 500 Instrument (Illumina). The quantification of transcript abundance was conducted using RSEM software (v1.2.22; supported by STAR aligner software (STAR 2.5.1b) and aligned to the hg38 human genome. Transcripts per million (TPM) values were then normalized among all samples using the upper quantile normalization method. Integration site features were calculated individually using each of the three RNA-Seq datasets from the three assayed participants; medians or averages of the three datasets were used for statistical analyses.

33A previously-described protocol with some modifications was used. Briefly, 20,000 total CD4 T cells were centrifuged at 1500 rpm for 10 min at 4°C in a precooled fixed-angle centrifuge. All supernatant was removed and a modified transposase mixture (including 25 μl of 2x TD buffer, 1.5 μl of TDE1, 0.5 μl of 1% digitonin, 16.5 μl of PBS, 6.5 μl of nuclease-free water) was added to the cells and incubated in a heat block at 37°C for 30 min. Transposed DNA was purified using a ChIP DNA Clean & Concentrator Kit (Zymo Research, #D5205) and eluted DNA fragments were used to amplify libraries. The libraries were quantified using an Agilent Bioanalyzer 2100 and the Q-Qubit™ dsDNA High Sensitivity Assay Kit (Invitrogen, #Q33230). All Fast-ATAC libraries were sequenced using paired-end, single-index sequencing on a NextSeq 500/550 instrument with v2.5 Kits (75 Cycles). The quality of reads was assessed using FastQC (https://www.bioinformatics.babraham.ac.uk). Low quality DNA end fragments and sequencing adapters were trimmed using Trimmomatic (http://www.usadellab.org). Sequencing reads were then aligned to the human reference genome hg38 using a short-read aligner (Bowtie2, http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) with the nondefault parameters “X2000”, “nonmixed” and “nondiscordant”. Reads from mitochondrial DNA were removed using Samtools (http://www.htslib.org). Peak calls were made using MACS2 with the callpeak command (https://pypi.python.org/pypi/MACS2), with a threshold for peak calling set to FDR-adjusted p<0.05. Integration site features were calculated individually using each of the three ATAC-Seq datasets from the three assayed participants; medians or averages of the three datasets were used for statistical analyses.

34We performed in situ Hi-C on sorted total CD4 T cells, central-memory and effector-memory CD4 T cells as previously described with minor modifications. Briefly, one million cells were crosslinked in 1% formaldehyde (Sigma, #F8775-25ML) for 10 min, followed by quenching with glycine (Sigma, #50046-50G). Cell nuclei were permeabilized with 0.4% sodium dodecyl sulphate (SDS) solution and chromatin was digested using 100 U of MboI (NEB, #R0147L). Overhangs generated by the restriction enzyme were filled using Klenow DNA polymerase (NEB, #M0210L) and a mix of dNTPs including biotin-14-dCTP (Invitrogen, #19518018). DNA fragments were then ligated using T4 DNA ligase (NEB, #M0202). After reversion of crosslinking, cellular proteins were digested with Proteinase K (NEB, #P8107), followed by DNA extraction using Phenol-Chloroform-Isoamyl alcohol. The extracted DNA was subjected to digestion of contaminating RNA with RNase A (ThermoFisher, #EN0531). Samples were sheared using a Covaris E220 instrument with a target size of 300-500 base pairs (140W peak incident power, each 67s, 10% duty, 200 cycles/burst) at the MGH NextSeq Core facility. Biotinylated fragments were isolated with Dynabeads MyOne Streptavidin C1 beads (Invitrogen, #65002). DNA linked to the Dynabeads was end-repaired using the End Prep Enzyme Mix and subsequently used for library construction with the NEBNext Ultra DNA Library Prep Kit for Illumina Sequencing (NEB, #E7370S). Final amplification of the libraries was done in 4 parallel reactions per sample according to the following program: 98 °C for 1 min, (98 °C for 10 s, 65 °C for 75 s, ramping 1.5 °C/s) repeated 12–20 times, 65 °C for 5 min, 4 °C. The indexed samples were processed individually and double size-selected using SPRISelect beads (Beckman Coulter, #B23318). Final Hi-C libraries were quantified using the Qubit dsDNA HS assay kit (Invitrogen, #Q33230) and a High Sensitivity D1000 kit (Agilent, #5067-5585) on a Tapestation 4200 (Agilent). The library was then sequenced on an Illumina NextSeq 500/550 (2x80 bp paired-end; NextSeq 500/550 High Output kit v2.5-150 cycles). HomerTools in the software suite HOMER (version 4.10.3) was used to truncate raw sequencing reads at the restriction enzyme cutting site, followed by aligning reads to the human reference genome (GRCh38) with Bowtie2 (version 2.3.4.3). The genome was tiled into bins with defined resolutions using the default method in HOMER to normalize the total read counts in each region. TADs were determined using default utilities from HOMER. Independently, we called significant inter-chromosomal and intrachromosomal interactions using FitHiC2 with Knight-Ruiz matrix balancing and FDR-adjusted p-values < 0.05 as a cutoff. Frequently interacting regions (FIREs) were called using FIREcaller (version 1.10). For data analysis, reads from autologous central-memory, effector-memory and total CD4 T cells were combined to generate a pooled dataset with overrepresentation of reads from memory cells, to account for the fact that most HIV-1-infected cells are included in the memory cell compartment. The numbers of total unique sequencing reads in these pooled datasets in the three assayed participants were 1,424,507,007; 1,685,072,487; and 2,400,898,418, respectively. Integration site features were calculated individually using each of the three Hi-C datasets; medians or averages of the three datasets were used for statistical analysis.

35Data are presented as circos plots, pie charts, bar charts, line graphs and box and whisker plots, showing the median, the 25% and 75% percentiles, and maximum/minimum values. Differences were tested for statistical significance using Mann-Whitney U tests (two-tailed), Fisher’s exact tests (two-tailed) or G-tests, as appropriate. Longitudinal correlations were evaluated using Kendall’s rank correlation coefficient (τ). p-values of <0.05 were considered significant. Analyses were performed using Prism (GraphPad Software, Inc.), Python (Python Software Foundation), and R (R Foundation for Statistical Computing).

36The Scikit-learn module in Python was used to develop a machine learning model to predict the transcriptional activity of HIV-1 integration sites. Briefly, the entire dataset of proviral species was randomly split into a training set with 80% of observations and a test set with 20% of observations, with each set containing the same proportion of observations from each participant and timepoint and the same proportion of RNA+ and RNA- proviruses. A logistic regression model was then trained on the training dataset with HIV-1 Long LTR RNA detection as the target variable. During training, stratified 10-fold cross-validation was used to determine the optimal data standardization and regularization strategies, and ultimately the model with the highest mean cross-validated area under the receiver operating characteristic curve (ROC AUC) was selected for evaluation with the test dataset.