HIV-seq reveals gene expression differences between HIV-transcribing cells from viremic and suppressed people with HIV

Domain: chemistry
Published: 2026-03-03
Editorial outcome: accept

Abstract
HIV-transcribing cells can perpetuate chronic inflammation in ART-suppressed people with HIV (PWH) and likely contribute to viral rebound after ART interruption. However, these cells are difficult to study using single-cell RNA-seq (scRNA-seq) due to their low frequency and low levels of HIV transcripts, which are usually not polyadenylated. By spiking in capture sequences targeting conserved regions of HIV during scRNA-seq – a new method we call “HIV-seq” – we detect double the mean number of HIV reads per cell from PWH. HIV RNA+ cells are enriched among T effector memory cells during both viremia and ART suppression but exhibit a cytotoxic signature during viremia only. In contrast, HIV-transcribing cells from ART-suppressed time points exhibit a distinct anti-inflammatory signature involving elevated TGF-β and diminished IFN signaling. These findings demonstrate that HIV-seq is a useful tool to better understand the mechanisms by which HIV-transcribing cells can persist during ART.

Introduction
Most people with HIV (PWH) experience rebound of HIV in plasma within several weeks after stopping antiretroviral therapy (ART), indicating the persistence of a rebound-competent viral reservoir that prevents HIV cure. The prevailing model has been that the rebound virus arises from a small fraction of HIV-infected cells that contain an infectious provirus in a latent state, where the infected cell does not constitutively produce virions but can be induced to do so after activation1–3. However, the rebound virus is often different from that which grows out ex vivo following stimulation in quantitative viral outgrowth assays (QVOA)4,5, suggesting that additional studies are needed to understand the reservoirs which can be reactivated in vivo.

A smaller body of research has focused on the cells actively transcribing HIV RNA in vivo, also known as the “active” reservoir. While HIV latency and expression are often viewed as a dichotomy (off/on), studies using multiple round QVOAs demonstrate varying degrees of inducibility ex vivo6, and blood and tissue cells from ART-suppressed individuals show a continuum of HIV-expressing cells in vivo, with variable degrees of progression through different blocks to HIV expression7–9.

Importantly, multiple lines of evidence suggest that HIV-infected cells which spontaneously express HIV RNA or protein in vivo may be just as important for pathogenesis and cure as the transcriptionally-silent reservoir. First, viral products expressed by active reservoir cells are likely to contribute to the immune activation and inflammation10,11 that are thought to underlie the sequelae of ART-treated infection, including organ damage and reduced life expectancy12–16. Second, the active reservoir seems better poised than the silent reservoir to immediately infect new cells after interruption of ART because it does not need to revert the mechanisms that prevent viral expression in the silent reservoir. Indeed, at least four studies have shown that different forms of cell-associated HIV RNA negatively correlate with time to rebound after ART treatment interruption (ATI)17–20. Moreover, one small study showed that in about half of people who interrupted ART, Pol sequences from the rebound virus matched those from cell-associated HIV RNA prior to ART interruption21. These findings support the hypothesis that some part of the active reservoir is rebound-competent. Finally, active reservoir cells, by expressing HIV RNA and/or protein, are likely more susceptible than quiescent latent cells to new immune-based therapies aimed at an HIV cure (e.g., therapeutic vaccines, TLR agonizts, adoptive immune therapies, and broadly neutralizing antibodies).

Despite all this evidence pointing toward the active reservoir as an important target for HIV cure, our understanding of these cells is still rudimentary, and it is unclear whether these cells exhibit similar phenotypic features as HIV RNA+ cells present during viremia. Cells transcribing or translating HIV genes from viremic individuals have been characterized by flow cytometry, CyTOF, and single-cell sequencing22–25 to a certain extent, but HIV-transcribing cells from ART-suppressed PWH have been harder to study with such single-cell technologies. For instance, conventional sequencing-based approaches identify active reservoir cells at such low numbers as to preclude meaningful analysis25. One alternative has been to activate cells from ART-suppressed individuals ex vivo to characterize reactivated HIV RNA+ cells26,27, but this approach characterizes ex vivo-stimulated and not spontaneously HIV-transcribing cells. Of note, sequencing-based approaches for detection of HIV RNA+ cells from virally-suppressed PWH are not only limited by the low throughput and high costs of droplet-based encapsulation technologies (e.g. 10X Genomics), but also the reliance on RNA capture through poly-dT. Because most HIV-infected cells from ART-suppressed individuals do not contain poly-adenylated HIV RNA due to blocks to transcriptional elongation and completion7, poly-dT-based methods of RNA capture will theoretically fail to recognize many HIV-transcribing cells.

To increase the ability to identify rare, HIV RNA+ cells from PWH, including in the context of ART suppression, we created a custom-modified, 10x Genomics-compatible, 5’ sequencing-based scRNA-seq workflow in which the poly-dT primer is supplemented with multiple HIV-specific reverse primers targeting different regions of the genome. This approach increases the capture of HIV-infected cells with non-polyadenylated HIV transcripts, including 5’ elongated as well as more distal HIV transcripts. In addition, we included DNA-barcoded antibodies (CITE-seq28) in our protocol to enable in-depth phenotyping of the HIV RNA+ cells for which we have transcriptome data. Compared to the standard 5’ sequencing approach, the inclusion of HIV-specific reverse primers allowed the detection of more HIV transcripts from PWH, and enabled for the first time a meaningful analysis of HIV RNA+ cells from ART-suppressed PWH. Using this advantage, we performed an in-depth analysis of the transcriptomes and phenotypes of HIV RNA+ cells from longitudinal samples of PWH during active viremia and after ART suppression.

Results
Development of HIV-seq as a method to increase capture and detection of HIV transcripts by single-cell sequencing
The standard 10X Genomics’ 5’ scRNA-seq workflow entails droplet encapsulation of individual cells, followed by capture and reverse transcription (RT) of polyadenylated transcripts using poly(dT) primers. This approach, theoretically, does not capture non-polyadenylated transcripts, although in practice some may be captured through non-specific binding. To more efficiently capture HIV transcripts, including those that are not poly-adenylated, we designed capture sequences targeting multiple conserved regions of the HIV genome (Fig. 1A): the R-U5-pre-gag region (for 5’ elongated transcripts), the pol gene (for mid-transcribed, unspliced transcripts), the second exon of tat-rev (for distal transcripts, including spliced), and two regions known to be enriched among intact proviruses: the packaging signal29 (Psi; for elongated, unspliced transcripts) and the HIV rev response element29 (RRE, for distal and unspliced/single-spliced transcripts). Our capture sequences match the consensus sequence of subtype B HIV-1, and are known to be efficient reverse primers for established ddPCR assays7. The HIV capture primers were spiked into the poly-dT primer mix and used in the RT following cell encapsulation (Supplementary Fig. 2). After RT, gene expression and CITE-seq libraries were prepared and sequenced. We aligned all sequences to the GRCh38 human genome, to which we had appended a subtype B HIV-1 consensus sequence we had generated from the Los Alamos database (see Methods). We named our overall pipeline “HIV-seq” due to its specific targeting of HIV transcripts for RT, library generation, and alignment.

To assess the utility of HIV-seq, we compared it head-to-head to the original 10X Genomics’ 5’ single cell RNA-seq pipeline (without HIV primers). To obtain sufficient numbers of HIV RNA+ cells for this comparative analysis, we selected samples from two viremic donors (PID1052 and PID8027) not on ART and processed equal numbers of cells using the two methods. Similar numbers of total CD4+ T cells were captured and analyzed with the original pipeline vs. with HIV-seq (PID1052: 4263 and 8772 cells, respectively; PID8027: 8022 and 9701 cells, respectively). We first assessed the extent to which HIV-seq may perturb the capture and sequencing of the host transcriptome. UMAP visualization revealed similar gene expression profiles in the absence vs. presence of HIV-seq (Fig. 1B). Differential gene expression (DEG) analysis in the absence vs. presence of HIV-seq revealed 4 DEGs in one donor and 6 DEGs in the other, with only a single gene (long non-coding RNA AL138963.4) overlapping (Supplementary Fig. 3, Supplementary Tables 3 and 4), suggesting minimal perturbations to human gene expression profiling upon HIV primer addition.

Next, we determined whether HIV-seq can lead to spurious detection of HIV RNA transcripts in HIV-uninfected cells. We first determined whether HIV-seq led to the identification of any HIV transcripts in non-permissive CD8+ T cells, and found that it did not (Supplementary Fig. 4A). We then performed a side-by-side comparison of HIV-seq vs. conventional 5’ scRNA-seq on PBMCs isolated from three HIV-uninfected donors. Similar numbers of total CD4+ T cells were captured and analyzed using the two methods (Donor 1: Conventional = 4668 cells, HIV-seq = 7260 cells; Donor 2: Conventional = 5417 cells, HIV-seq = 5998 cells; Donor 3: Conventional = 12,327 cells, HIV-seq = 8753 cells). Not a single HIV transcript was identified, irrespective of HIV capture primers (Supplementary Fig. 4B). Given the lack of background signal, this experiment also established that one HIV transcript could be used as the threshold for defining HIV RNA+ cells using HIV-seq.

Quantitative analyses showed that HIV-seq identified a higher percentage of HIV RNA+ cells from both donors than did traditional 5’ scRNA-seq (Fig. 1C). For PID1052, HIV-seq increased the identification of HIV RNA+ cells from 0.047% (2 HIV RNA+ cells in 4263 CD4+ T cells) to 0.068% (6 HIV RNA+ cells in 8772 CD4+ T cells), corresponding to a 44.7% increase. For PID8027, HIV-seq increased the capture of HIV RNA+ cells from 0.77% (62 HIV RNA+ cells in 8022 CD4+ T cells) to 0.97% (93 HIV RNA+ cells in 9701 CD4+ T cells), corresponding to a 26% increase. While the small number of study participants did not support formal statistical comparison of HIV RNA+ cell frequencies, HIV-seq also detected a significantly higher number of HIV transcripts per infected cell (mean of 48 vs. 25 HIV reads/cell; P = 0.0012; Fig. 1D).

To assess whether the additional HIV transcripts captured by HIV-seq preferentially mapped to certain regions of the HIV genome, we compared the distribution of the HIV reads from the two viremic individuals in the absence vs. presence of HIV-capture primers. HIV-seq increased the detection of HIV transcripts from across the entire proviral genome, without changing the representation of the HIV regions detected relative to traditional 5’ scRNA-seq (Fig. 1E). Under both conditions, most HIV transcripts aligned to the pol and early gag regions of HIV. Application of HIV-seq in the context of ART suppression also resulted in detection of transcripts primarily mapping to the gag and pol regions (Supplementary Fig. 5). These results demonstrate that HIV-seq detects more HIV transcripts than conventional 5’ scRNA-seq and suggest that HIV-seq enables the detection of more HIV RNA+ cells.

HIV RNA+ cells from viremic PWH are preferentially cytotoxic Tem cells and exhibit an intracellular state promoting HIV replication
Leveraging the ability of HIV-seq to increase the numbers of HIV RNA+ transcripts and, to some extent, the number of HIV RNA+ cells we can analyze by scRNA-seq, we implemented it on cells from 4 viremic, ART-naïve PWH. We identified 8, 869, 155, and 40 HIV RNA+ cells out of 13,035, 35,849, 17,723, and 19,060 CD4+ T cells from PID1052, PID8026, PID8027, and PID145, respectively. In total, 1072 HIV RNA+ CD4+ T cells were identified out of the 85,667 CD4+ T cells analyzed; infected cell frequencies ranged from 0.061% to 2.42% of the CD4+ T cell population (Fig. 2A). UMAP visualization of the transcriptomes of the HIV RNA+ cells revealed heterogeneity, in that infected cells were distributed in multiple regions of the UMAP space, yet enrichment was observed in some regions, suggesting non-random distribution of infected cells among CD4+ T cell subsets (Fig. 2B).

Considerable variability in the degree of HIV transcription was observed among infected cells, with HIV transcript levels ranging from 1 to 1063 HIV reads per cell. When we separated the HIV RNA+ cells into those with low numbers of HIV reads (HIV low: 1 to 50 HIV reads) and those with high numbers of HIV reads (HIV high: >50 HIV reads), we found that participants with lower numbers of HIV RNA+ cells (PID1052 and PID0145) only harbored HIV low cells, while those with higher numbers (PID8026 and PID8027) had both populations (Supplementary Fig. 6A). Of note, HIV-seq did not affect whether HIV high could be detected (Supplementary Fig. 6B). The distribution of HIV high and HIV low cells across the UMAP was similar (Supplementary Fig. 6C) and the only significantly differentially expressed transcript/protein between these populations – apart from HIV RNAs – was the CD4 protein, which was decreased among the HIV high cells (Supplementary Fig. 6D). This finding likely reflects higher expression of Nef, which downregulates cell-surface CD430, in the HIV high cells. Because the HIV high and HIV low cells exhibited overall similar gene expression profiles, all remaining analyses combined these two populations together.

We then assessed whether the HIV RNA+ cells were enriched in specific cellular subsets. We first assessed the distribution of classical CD4+ T cell subsets (Supplementary Fig. 7) among uninfected and HIV RNA+ CD4+ T cells. Relative to their uninfected CD4+ T cell counterparts, HIV RNA+ cells were under-represented among naïve cells and over-represented among memory cells (Fig. 2C), as expected31. Within memory T cells, HIV RNA+ cells were under-represented among central memory (Tcm) cells and those of the CCR7+ CD27- phenotype, and were over-represented among effector memory (Tem) cells (Fig. 2C), consistent with prior reports25,32. By contrast, transitional memory (Ttm), regulatory T cells (Treg), and T follicular helper (Tfh) cells were equally represented among uninfected and HIV RNA+ CD4+ T cells (Fig. 2C).

Next, we compared differentially expressed proteins (DEPs) and transcripts between HIV RNA+ and HIV RNA- cells from the viremic PWH. At the protein level, HIV RNA+ cells expressed lower levels of the naïve T cell marker CD45RA and higher levels of the memory T cell marker CD45RO compared to HIV RNA- cells across all donors (Supplementary Fig. 8), consistent with the notion that HIV RNA+ cells are enriched in memory CD4+ T cells (Fig. 1C). Other differentially expressed proteins were only observed in a subset of donors; these included CD279 (PD1; upregulated in HIV RNA+ cells in 2/4 donors), CD27 (downregulated in HIV RNA+ cells in 3/4 donors), and CD49d (upregulated in 3/4 donors) (Supplementary Fig. 8).

At the RNA level, almost 300 genes were differentially expressed between HIV RNA+ and HIV RNA- cells from the viremic PWH (Fig. 2D, Supplementary Table 5). Consistent with disenrichment of Tcm among HIV RNA+ cells (Fig. 2C), transcript levels of the CCR7, SELL (encoding the protein CD62L), and CD27 – markers of Tcm cells33,34 – were downregulated in the HIV RNA+ cells (Fig. 2D). HIV RNA+ cells also had low transcript levels of the alarmins S100A8 and S100A9, which encode for S100 calcium binding protein A8 and A9, respectively. These proteins are released in response to environmental triggers and cellular damage35. Antiviral factors, including SERPINA136,37 and APOBEC3A38, were also decreased among HIV RNA+ cells, suggesting that the intracellular state of HIV RNA+ cells in viremic individuals may favor HIV replication. In line with this finding, HIV RNA+ cells also exhibited downregulation of CST3, which encodes cystatin C, a cysteine protease inhibitor that interacts with the HIV proteins gp160, gp120, p31 and p24, and inhibits HIV protease function39. With regard to transcripts upregulated among HIV RNA+ cells, HIV transcripts were the top hit, as expected. In addition, HIV RNA+ cells expressed higher levels of the activation marker CD70. Interestingly, CD4+ CD70+ cells are increased in PWH with high levels of viremia and associate with immune activation40, suggesting that this subset of infected cells may contribute to disease progression during untreated infection. HIV RNA+ cells also expressed elevated levels of CXCR6, which encodes a chemokine receptor that is an alternative co-receptor for HIV41.

Pathway analysis of the DEGs (Fig. 2E) revealed upregulation of the NFAT pathway (MAF, CLTA4) in HIV RNA+ cells, consistent with NFAT as a driver of HIV transcription42,43. Similarly, there was an upregulation of the PKC pathway, known for its involvement in HIV gene expression and latency reversal44. Finally, and consistent with the DEG analysis, chemokine signaling pathways – featuring genes, such as CXCR6, the CXCR6 ligand CXCL13, and PLCB1 – were also elevated among HIV RNA+ cells.

While DEG analysis identified both known and novel shared features among all HIV RNA+ cells, it was clear that the HIV RNA+ cells were heterogeneous (Fig. 2B). We therefore implemented clustering to assess for transcriptomic or phenotypic features that may not be shared by the entire population of infected cells. Louvain clustering identified thirteen clusters of CD4+ T cells (Fig. 2F). The classical CD4+ T cell subsets (Supplementary Fig. 7) did not define the clusters (Supplementary Fig. 9A), and in fact all the classical memory subsets were represented among the 13 clusters, albeit in different proportions (Supplementary Fig. 9B). The 13 clusters could be differentiated by lineage markers, effector and cytotoxic profiles, and activation/proliferative states (Supplementary Fig. S10). HIV RNA+ cells were differentially distributed among the 13 clusters as compared to their uninfected counterparts, with highest enrichment of HIV RNA+ cells in cluster 6; 37% of all HIV RNA+ cells were located in this cluster, as compared to only 17.3% of total CD4+ T cells and 17.1% of HIV RNA- cells (Fig. 2F and Supplementary Table 6). Cluster 6 comprised mostly Tem cells (Supplementary Fig. S9B), and was characterized by high expression of cytotoxic and cytolytic genes, including GZMA, GZMB, GZMH, GZMM, PRF1, GNLY, and NKG7, as well as the Th1-defining factors IFNG and TBX21 (Fig. 2G and Supplementary Fig. 10). This finding suggests that cells from cluster 6 are Tem cells of a cytotoxic Th1 phenotype and is consistent with prior reports that HIV RNA+ cells from viremic PWH exhibit Th1 cytotoxic signatures25,32.

Overall, our data demonstrate that HIV RNA+ cells from viremic PWH are heterogeneous but exhibit shared features, including being more likely to be Tem, and notably, displaying a previously undescribed state more conducive to viral replication, with low expression of restriction factors and increased activation of cellular pathways promoting HIV gene expression. Additionally, relative to their uninfected counterparts, they preferentially exhibit a cytotoxic signature characterized by higher expression of granzymes, perforin, and granulysin and a Th1 signature.

HIV RNA+ cells from suppressed PWH are preferentially Tem cells but do not exhibit a cytotoxic signature
Three of the individuals we had analyzed in the context of active viremia (PID1052, PID8026, and PID8027) had specimens collected after >24 weeks of ART suppression. Therefore, we next implemented a similar HIV-seq analysis pipeline to characterize these participant-matched HIV RNA+ cells in the context of ART suppression. We identified 14, 7, and 4 HIV+ cells out of 15,468, 42,618, and 17,072 CD4+ T cells for PID1052, PID8026 and PID8027, respectively. A total of 25 HIV RNA+ cells were detected out of the 75,158 CD4+ T cells analyzed, corresponding to an HIV RNA+ cell frequency of 0.016% to 0.091% (Fig. 3A). As for the viremic specimens, HIV RNA+ cells were broadly distributed among CD4+ T cells, demonstrating heterogeneity (Fig. 3B), and were enriched among Tem cells (Fig. 3C). In fact, in two individuals (PID8026 and PID8027), the entire HIV RNA+ cell population was exclusively found within Tem cells. When comparing HIV RNA+ and HIV RNA- cells, we found no DEPs that were consistent across study participants (Supplementary Fig. 11). Likewise, few DEGs were observed between the HIV RNA+ and HIV RNA- CD4+ T cells, with the notable exception of a handful of downregulated genes (Fig. 3D and Supplementary Table 7). These included the Tcm marker SELL, the alarmin S100A9, and the cysteine protease inhibitor CST3, all of which were also preferentially downregulated among HIV RNA+ cells during active viremia (Fig. 2).

Louvain clustering analysis revealed that unlike HIV RNA+ cells collected during viremia, those collected during suppression were not enriched among cluster 6 (cytotoxic Tem). Instead, the cluster distribution of the infected cells under ART revealed that HIV RNA+ cells were primarily distributed among cluster 1, with 52% of all HIV RNA+ cells being in cluster 1 as compared to only 28.1% in both total CD4+ T cells and HIV RNA- cells (Fig. 3E and Supplementary Table 6). Cluster 1 comprises CD4+ T cells expressing high levels of the IL7R transcript (Supplementary Fig. 10), which encodes for a receptor that is characteristic of long-lived cells and is important for homeostatic proliferation. Interestingly, cluster 1 also includes cells with Th2 features as defined by high GATA3 expression (Supplementary Fig. 10), although IL4 and IL5 – cytokines characteristics of Th2 cells – were not detected.

Our observation that cells in both clusters 1 and 3 express high levels of IL7R support the notion that the active HIV reservoir persists in long-lived/stem-like cells with proliferative potential. Indeed, IL7R expression – at both the transcript and protein levels – was spatially localized to regions of the UMAP occupied by clusters 1 and 3 (Fig. 3F). Furthermore, when we directly compared the HIV RNA+ cells to HIV RNA- cells, we found that IL7R transcript and IL7R protein were both more highly expressed on the HIV RNA+ cells (P < 0.05; Fig. 3G). That HIV RNA+ cells during ART persist in stem-like cells is further supported by our observation that a substantial proportion (44%) of the HIV RNA+ cells express BCL-2, a pro-survival, anti-apoptotic gene implicated in driving long-term persistence of HIV46–49 (Fig. 3H). Taken together, these findings suggest that in blood, most HIV RNA+ cells during ART suppression reside not within cytotoxic CD4+ T cells – as was observed during viremia – but rather within a long-lived population of memory CD4+ T cells.

CD4+ T cells exhibit stronger antiviral immunity and upregulate the integrated stress response pathway during viremia as compared to after ART suppression
Leveraging the fact that HIV-seq was performed on paired specimens from before vs. after ART suppression, we then compared the transcriptomes of total CD4+ T cells across these two conditions. Total CD4+ T cells coming from viremic vs. suppressed specimens were transcriptionally divergent, as reflected by distinct UMAP localizations (Fig. 4A, Supplementary Fig. 12), even though their distribution among classical CD4+ T cell subsets was similar (Fig. 4B). Although DEGs were identified (Fig. 4C, Supplementary Table 8), none retained statistical significance following correction for multiple comparisons. As an exploratory analysis, however, we assessed the DEGs with the lowest raw p-values to gain insights into potential cellular pathways distinguishing CD4+ T cells during active viremia from those following ART suppression. This analysis revealed that during viremia, CD4+ T cells preferentially express higher levels of cytotoxic and proinflammatory genes (SP100, IL6), and genes related to interferon α, β and γ signaling (IFI27, IFI44L, IFIT3, IFI6, IFI44, ISG15, IFITM10, IFI30, IFI27L1, MX1, IFIH1, IFNLR1, IFIT5, NEAT1, TYMP). Some of these (IFI44L, ISG15, NEAT1, TYMP) have previously been reported to be higher among CD4+ T cells during active viremia as compared to ART treatment50. This upregulation may reflect an antiviral host response triggered by high levels of viral transcripts and proteins during viremia. Other highly expressed genes during viremia included cytotoxic genes (GZMB, GZMK) whose expression can be induced by cytokines, such as IL2 and IL1551, and during inflammation and viral infections52.

Also of interest was that during viremia, CD4+ T cells increased expression of EIF2AK2 (Fig. 4C, Supplementary Table 8), a key gene involved in integrated stress response (ISR), a pathway previously reported to be induced during acute HIV infection53,54. During ISR, EIF2AK2 expression is induced by type I IFNs, which then upon binding to viral dsRNA can initiate a cascade of events culminating in diminished protein translation. Indeed, when we performed DEG analysis of all donors combined using an approach previously implemented to identify DEGs between total CD4+ T cells during active HIV viremia vs. after ART suppression25,50 (see Methods), we observed downregulation of numerous ribosomal transcripts (RPL and RPS transcripts) during viremia (Supplementary Fig. 13 and Supplementary Table 9). The diminished ribosomal transcript levels were accompanied by diminished expression of EEF2, whose downregulation has been shown to reduce active ribosomes54 and shutdown mRNA translation, resulting in overall diminished viral protein synthesis55. Hence, the ISR pathway, which serves to coordinate cellular responses to various stressors by regulating protein synthesis and gene expression and is a host response to limit viral spread, appears to be a characteristic feature of acute untreated HIV infection that is turned off upon ART suppression.

Relative to HIV RNA+ cells during viremia, those during ART suppression upregulate TGFβ signaling pathways and exhibit diminished activation of host responses
We then compared the HIV RNA+ CD4+ T cells from the paired viremic and ART-suppressed timepoints. Classical CD4+ T cell subset distributions of HIV RNA+ cells were similar between the two timepoints, with preferential distribution among Tem cells for both (Fig. 4D). DEP analysis revealed lower CD4 expression on HIV RNA+ cells during viremia as compared to during ART suppression (Supplementary Fig. 14), perhaps due to elevated expression of Nef, an HIV accessory gene that potently downregulates cell-surface CD4 expression in productively-infected cells56. Multiple DEGs were identified between the timepoints, but none retained statistical significance after correction for multiple comparisons. As an exploratory analysis, we assessed the identities of the DEGs based on the raw p-values (Fig. 4E and Supplementary Table 10). This analysis revealed upregulation of multiple ISGs and proinflammatory genes in HIV RNA+ cells during viremia compared to suppression, including IFI27, MX1, ISG15, DUSP2, and SP100. These genes were also upregulated among total CD4+ T cells during viremia (Fig. 4C and Supplementary Table 8). IFI27 in particular was increased by more than 10-fold in HIV RNA+ CD4+ T cells during viremia as compared to those during ART. Given that HIV-1 Vpr and Tat can directly induce IFI27 production57–59, and that expression of these viral proteins may be diminished in the context of ART suppression due to multiple blocks to HIV transcription7–9, it is possible that elevated expression of Vpr and Tat is driving IFI27 expression prior to ART initiation. Interestingly, IFI27 levels correlate with inflammation and disease progression during both HIV-1 and HIV-2 infection60,61, and they may contribute to pathogenesis60,61.

DEGs associated with HIV transcription were also upregulated in HIV RNA+ cells during viremia as compared to after ART suppression. HIV RNA+ cells during viremia expressed higher levels of SRRM1, a modulator of HIV-1 splicing that is involved in the regulation of Tat and Nef expression62, and lower levels of two genes implicated in silencing HIV transcription and translation: RBL2, which recruits and targets histone methyltransferases, leading to epigenetic transcriptional repression; and RPS10, which binds HIV Nef to form a complex that decreases viral protein synthesis63. Together, these findings suggest that during viremia, HIV RNA+ cells exhibit a transcriptional profile that favors the production of more HIV transcripts. This finding accords with our having observed elevated HIV transcript levels among HIV RNA+ cells during viremia as compared to during suppression in two out of our three participants (11- and 14-fold increase, respectively, for PID8026 and PID8027).

Pathway analysis of the DEGs in HIV RNA+ cells before versus after ART suppression supported the notion that HIV RNA+ cells during viremia are preferentially in an activated, antiviral state, characterized by upregulation of multiple interferon pathways (IFN I and IFN II) (Fig. 4F, top). Conversely, HIV RNA+ cells during ART suppression preferentially upregulated TGF-β signaling (Fig. 4F, bottom). This finding was driven by upregulation of RBL2 and ITGB1 (Supplementary Table 10), genes associated with the TGF-β signaling pathway. ITGB1- and RBL2-associated TGF-β activation has been implicated in tumor suppression and cancer growth arrest64,65. Although ITGB1 and RBL2 have not been directly implicated in HIV infection, TGF-β signaling has been recently implicated in promoting HIV latency66,67. Our data suggest that this TGF-β-associated signature is a phenomenon that only emerges in the context of ART suppression, as HIV RNA+ cells during viremia do not exhibit such a signature. Overall, these results indicate that the transcriptional profiles of HIV RNA+ cells differ depending on whether or not ART is present, and that the features of HIV RNA+ cells during viremia cannot be assumed to be the same as those during ART suppression.

Methods
Ethics statement
The study was approved by the Committee on Human Research (CHR), the Institutional Review Board for the University of California, San Francisco (approval #11–07551 and #10-01561). All study participants provided written informed consent.

Study Population
The study participants were 4 people living with HIV (median age = 38.5; median CD4 count = 490 cells/mm3; Table 1). Participants were selected based on availability of cryopreserved PBMC from before the start of ART (n = 4), with a preference for those with an additional longitudinal sample of cryopreserved PBMC obtained after ART suppression (n = 3). Paired and longitudinal, archived PBMC samples were obtained from the UCSF Treat Acute Study and the SCOPE cohort. Samples were collected prior to ART initiation (Week 0) and following ART suppression (Week 24: PID8026; Week 45: PID8027; and Week 70: PID1052). A total of two aliquots of 107 cells each were obtained for each time point (viremic and suppressed), with one aliquot preserved for HIV DNA/RNA measurements. An additional participant (PID0145), for whom only the viremic (Week 0) time point was available, was recruited from the San Francisco VA Medical Center. Three additional PBMC samples from HIV-uninfected donors were used as negative controls, to establish the extent to which HIV-seq results in false-positive identification of HIV reads and to set up the HIV transcript threshold for identification of HIV-infected cells.

In this pilot study where it was only feasible to study a small number of participants, and given the difficulty in finding paired longitudinal samples from before and after ART, we accepted the first available samples, which were provided to us in de-identified form. Therefore, sex/gender was not considered in the study design. Data on sex/gender was obtained after completion of data analysis, and was self-reported. With only three male participants and one female participant, the numbers were deemed insufficient for a post-hoc analysis.

HIV DNA measurements
For each participant, total HIV DNA (R-U5-pre-gag region) levels were measured to ascertain whether their reservoir was sufficiently high for detection with limited cell inputs (up to 10,000 cells for each well of the 10X Genomics scRNA-seq platform). An aliquot of 107 cryopreserved cells was tested from each time point (viremic and suppressed) for each individual. Each aliquot was thawed, cryopreservation medium was removed, and total RNA and DNA were extracted using TRI Reagent (Molecular Research Center, Inc., Cincinnati, OH) as per manufacturer’s instructions, with the following modifications: polyacryl carrier (Molecular Research Center, Inc., Cincinnati, OH) was added to TRI reagent prior to lysis, RNA was resuspended in RNase free-water, DNA was extracted using back extraction buffer (4 M guanidinethiocyanate, 50 mM sodium citrate, 1 M Tris), polyacryl carrier was added to the aqueous phase containing the DNA, and DNA was resuspended in QIAGEN buffer EB8. Cellular DNA was fragmented using a QIAshredder column (QIAGEN) to facilitate incorporation into droplets.

HIV DNA and copies of the housekeeping gene TERT (telomere reverse transcriptase) were measured in duplicate by droplet digital (dd)PCR (Bio-Rad QX200)7. Each reaction consisted of a 20 µL solution. The ddPCR for TERT contained 10 µL of ddPCR Probe Supermix (no dUTP) [Bio-Rad], 1 uL of 20x TaqMan™ Copy Number Reference Assay, human, TERT (Thermo Fisher Scientific), and 50 ng of cellular DNA. The ddPCR for HIV DNA contained 10 µL of ddPCR Probe Supermix (no dUTP) [Bio Rad], 900 nM of each primer (F: gcctcaataaagcttgccttga; R: gggcgccactgctagaga; both from Thermo Fisher Scientific), 250 nM of probe (5’ FAM-ccagagtcacacaacagacgggcaca-3’ nonfluorescent quencher/Minor Groove Binder [MGB]; Thermo Fisher Scientific), and 500 ng of cellular DNA. Droplets were amplified with the following cycling conditions: 10 min at 95°C, 45 cycles of 30 s at 95°C and 59°C for 60 s, and a final droplet cure step of 10 min at 98°C. Droplets were read and analyzed using the QuantaSoft software in the “Absolute” quantification mode. A threshold of three HIV DNA copies/10,000 cells at the suppressed time point was set as a threshold for study inclusion, as this infection frequency was predicted to yield a reasonable likelihood of detecting at least one HIV RNA+ cell using the 10X Genomics scRNA-seq platform. All samples subjected to single-cell sequencing in this study met this threshold.

scRNA-seq sample preparation
PBMC. A total of 107 PBMCs from each donor were thawed at 37°C and washed once with warm media (RPMI [Corning Inc., Corning, NY] supplemented with 10% FBS [VWR, Radnor, PA]) and then resuspended in FACS buffer (RPMI, supplemented with 2% FBS and 2 mM EDTA [Thermo Fisher Scientific]) prior to counting. After setting aside 105 cells from each sample for downstream PBMC spike-in, CD4+ T cells were purified from the remaining cells by negative selection using the EasySep Human CD4+ T Cell Enrichment Kit (StemCell). The PBMCs from the same donor that were set aside were then spiked back in a ratio of 5:100 in order have non-CD4+ T cells to help establish gates for defining HIV RNA+ cells.

TotalSeq-C Antibody Staining. TotalSeq-C pooled antibody mix (Biolegend) was prepared so as to contain 0.4 μg of each antibody in a final volume of 100 μL/sample, in PBS (Ca++ and Mg+ free, Corning Inc.) containing 3% FBS, as per the devised panel (Supplementary Table 1). A million cells were pelleted, resuspended in RPMI supplemented with 3% FBS, and incubated at 4°C for 10 min with Fc block (Miltenyi) at a 1:10 dilution. 100 μL/sample of TotalSeq pooled antibody mix was then added to the cells, which were incubated for 30 min at 4°C. Cells were then washed three times in RPMI containing 0.04% BSA, strained using a 40 μm cell filter (BD Falcon), and resuspended at a concentration of 1000 cells/μL. The TotalSeq-C panel was used to facilitate manual gating: CD45 to identify immune cells; CD19, CD20, CD14, and CD8 to exclude B cells, myeloid cells, and CD8+ T cells; CD3 and CD4 to identify CD4+ T cells; CD45RA and CD45RO to distinguish between naïve and memory CD4+ T cells; and CD197 (CCR7), CD27, CD62L, CD127, CD279 (PD1), and CD49d as additional subset/phenotype markers.

Custom HIV primers
Custom HIV-seq primers (obtained from Thermo Fisher Scientific) were designed by appending HIV-specific capture sequences (Supplementary Table 2) to a non-poly(dT) PCR handle (Fig. 1). The concentration of 10X Poly(dT) oligo (poly-dT RT Primer PN 2000007) is estimated to be 2.7 μM (100 pmol in 36.8 μl of RT mix) in each 10X RT reaction93. Initial experiments spiking in different concentrations of pre-gag primers revealed that 0.67 μM of primer resulted in a significantly lower total cDNA yield as compared to the 0.33 μM condition and the condition where no HIV primers were added (Supplementary Fig. 1). Hence, we implemented 0.33 μM of the primers moving forward. The primer pool (containing 84 nM of each primer) was added to the 10X Genomics RT reaction mix containing RT Reagent B, Poly-dT RT Primer, Reducing Agent B, and RT Enzyme C (HIV-Seq RT Mix), which was used for RT of encapsulated cells as described below.

RT, library preparation and sequencing
Chromium Next GEM Single Cell 5’ Reagent Kits v2 (Dual Index) with Feature Barcode technology for Cell Surface Protein & Immune Receptor Mapping (10X Genomics, PN 1000263) and the Chromium™ Controller (10X Genomics, PN 120223) were used for gene expression library and cell-surface protein library generation. Briefly, 20 μl of a TotalSeq C-stained single-cell suspension (at ~1000 cells/µl, which corresponds to ~20,000 cells) was mixed with the HIV-Seq RT Mix, barcoded using Single Cell VDJ 5’ Gel Beads, and partitioned with oil onto a Chromium Next GEM Chip K, following the 10X Genomics protocol for a targeted recovery of 10,000 cells. The chip was then loaded onto a Chromium Controller for single-cell GEM generation and barcoding. RT reactions were performed according to the Chromium Single Cell 5’ Reagent Kits v2 (Dual Index) User Guide (10x Genomics CG000330, Rev A). Sequencing libraries were constructed with 13 cycles of PCR during cDNA amplification and 14 cycles of Sample Index PCR. Gene expression and cell surface protein libraries were pooled and sequenced on a NovaSeq 6000 lane (S4 flow cell) and sequenced at a minimum of 50,000 reads/cell.

Data processing, statistics, and analysis
Sequencing libraries from n=4 independent experiments were analyzed. A custom Human/HIV consensus subtype B reference sequence was purpose-built for this study by downloading a Group M alignment from the Los Alamos National Laboratory’s HIV sequence database (2018 version), filtering for subtype B sequences, realigning using MAFFT94 with default settings, and generating a majority consensus sequence using Geneious (http://www.geneious.com/). The consensus HIV sequence was appended to the human reference genome (GRCh38-2020) and annotated as a single exon. Alignment of scRNA-seq reads, collapsing of reads to unique molecular identifier (UMI) counts, and cell calling was performed using CellRanger 6.0.2 (10X Genomics). Filtered count matrices of features generated using the CellRanger count function were then subjected to multiple cleanup steps. Cells with high mitochondrial gene expression (mtDNA% > 15%) and low features (<500) were removed, along with doublet cells (identified using the DoubletFinder package95), in R. All samples were normalized for their transcriptome library depth and batch-corrected between donors using the Seurat96 NormalizeData and integration function from Seurat package, in R. Harmony batch correction was applied for visualization where indicated. Downstream analysis was performed in SeqGeq (mostly for visualization; FlowJo, LLC) and Seurat v4.3.0. In all analyses, genes are depicted in italics and proteins in bold.

Clustering. Graph-based clustering was performed using the Louvain algorithm implementation97 in the FindClusters Seurat function. The optimal clustering resolution parameters were determined using Random Forests98 and a silhouette score-based assessment of clustering validity and subject-wise cross-validation99. This approach entails evaluating the clustering of cells in a given biological sample based on predictions derived after holding out cells from this sample, and clustering the cells from all the remaining samples, and ensures reproducible clusters across the biological samples. Clusters that appeared biologically under-clustered were subsequently re-evaluated using the same optimal resolution determination procedure and sub-clustered. T cell Receptor Alpha Variable (TRAV) and T cell Receptor Beta Variable (TRBV) genes were removed from the variable features used for clustering as they were driving the clustering in a donor-specific manner, as expected with randomly-generated VDJ sequences. HIV was also removed from the variable features used for clustering. Thirteen distinct biologically relevant clusters were identified, annotated manually, and used for further analyses.

Manual gating. Manual gating was conducted using SeqGeq software to delineate classic CD4+ T cell subsets and to identify HIV-expressing cells. Both gene and protein expressions were used for gating.

Statistical Analysis of CD4+ T Cell Subsets and HIV RNA+ Cell Distribution. For establishing associations between CD4+ T-cell subset proportions among the HIV RNA- and HIV RNA+ cells, and HIV RNA+ cell proportions among the different clusters, we used a generalized linear mixed model (GLMM) with a binomial probability distribution implemented in the lme4100 package in R. In the model, CD4 subset or cluster membership was treated as the outcome being studied. This outcome was represented by comparing the number of cells within the subset/cluster to the number of cells outside of it. The change in subset/cluster membership between HIV RNA status and timing of measurements (viremic vs. suppressed) was estimated as a log odds ratio, defined as the change in the log odds of subset/cluster membership. This change was estimated with the emmeans101 R package using the GLMM model fit.

Differential expressed protein (DEP) analysis. Raw ADT counts were normalized to 10,000 reads per cell to obtain relative expression values for each marker per cell. Marker expression was then visualized using violin plots, which display data density, median, and interquartile range, to compare the two populations of interest (as indicated in the figure legends) for each donor. Statistical significance was assessed using the Mann–Whitney U test for each marker, followed by False Discovery Rate (FDR) correction of p-values to adjust for multiple comparisons.

Differential expressed gene (DEG) analysis. Pseudobulk DEG analysis was performed to identify genes or proteins that were significantly upregulated or downregulated between different cell populations or experimental conditions. DEG analysis at the RNA level was carried out using the R package muscat102, and counts were summed across clusters and samples. The pbDS function to estimate associations with disease was run with the edgeR103 method, with minimum cells set to 3 and filter set to gene, and donor was included as a confounder in the model. In the exploratory analysis comparing total CD4+ T cells during viremia vs. suppression, all individuals within a group (viremic or ART-suppressed) were combined, and then gene expression levels were compared across individual cells within each group, followed by implementation of the Wilcoxon rank-sum test to assess for DEGs and proteins, similar to analytical approaches recently described25,50. Of note, this analytical approach should be considered exploratory as it does not fully account for correlations between cells from the same subject or paired study designs104,105. Results visualized as volcano plots were plotted using EnhancedVolcano106. Select genes and proteins of interest among parameters that passed the indicated adjusted p-value (or raw p-value where indicated) thresholds of 0.05 and a log2 fold-change greater than 0.25 are highlighted in the volcano plots. All genes are depicted in italics and proteins in bold. T cell receptor (TCR) genes (TRAV and TRBV) were deliberately omitted from the analysis, since private TCR sequences drove donor-specific effects, and the primary focus of our study was to identify gene expression profiles shared across donors.

Pathway Enrichment Analysis. The differential gene expression analysis results filtered by adjusted p-value < 0.05 and log2 fold-change > 0.25 were subjected to over-representation analysis (ORA) using the Enrichr107–109 web-based tool. Pathway analysis was performed against two curated databases: BioPlanet 2019 and MSigDB Hallmark 2020. This analysis identified significantly enriched pathways and functional categories along with statistical metrics and p-values.

Distribution of HIV reads
An additional analysis assessing distribution of reads across the HIV genome was performed to determine where HIV transcripts most frequently align. We selected PID8027 and PID1052, for whom we had data on the absence versus presence of HIV-capture primers. Gene expression reads from the week 0 (viremic) timepoint were aligned to a custom, annotated HXB2 reference sequence, which defined discrete HIV coding regions (5’ LTR U3, 5’ LTR R, 5’ LTR U5, start of gag, gag P1, pol, vif, vpr, rev, vpu, env, tat, nef, 3’ LTR U3, 3’ LTR R) using CellRanger (10X Genomics, version 6.0.2). The distribution of HIV reads aligning to these specific regions was then filtered, quantitated using an R script, and compared between samples subjected to the conventional 10X Genomics’ 5’ scRNA-seq pipeline vs. HIV-Seq.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this paper, Frouard et al. developed a new single-cell sequencing method, HIV-seq, to perform targeted enrichment of HIV transcripts and improve the efficiency for detecting the rare HIV RNA+ cells from PWH. The authors incorporated this method with single-cell CITE-seq and profiled clinical samples from PWH before and after ART treatment. By comparing the phenotypes during viremia and ART, the author described the signatures of HIV-infected cells during viremia and after viral suppression. The major strength of this study is the adding spike-in primers before cells are loaded into the 10x Chromium machine for targeted enrichment increased the capture of HIV RNA+ cells which increased HIV capture rate. The iterative adjustment of spike-in primers is a good endeavor optimizing this assay. The major weakness is the lack of demonstration of sensitivity, specificity, and how these spike-in primers biased cellular transcriptome. The lack of threshold of setting HIV RNA+ reads make the data not trustable, and thus an uninfected donor control is needed to set the appropriate threshold and to avoid false positive calling of HIV+ cells. The cell cluster annotation and differential gene expression analysis are not rigorous and should be improved. Having p <= 0.05, not p <0.05, is a fundamental statistical flaw. If well addressed, this study can bring innovation and value for the field.\nMajor comments:\n1. For any new single-cell assay, a sensitivity and specificity assay, and how these spike-in primers biased cellular transcriptome must be established first. Please primary CD4+ T cells infected with HIV and sorted/enriched (to obtain 100% HIV infection) versus primary CD4+ T cells with no infection, to perform HIV-seq versus regular 10x 5’ CITEseq [4 samples: (a) sorted 100% infected CD4 with HIV-seq, (b) sorted 100% infected CD4 with regular 5’ CITEseq, (c) uninfected CD4 with HIV-seq, (d) uninfected CD4 with regular 5’ CITEseq]. This benchmarking experiment is needed to address: (a) How should the authors set the threshold of HIV detection in uninfected cells, such as using HIV read count or UMI count as a threshold? (b) what’s the sensitivity detecting HIV RNA+ cells increase by using HIV-seq?\n2. Impact of spike in on transcriptome: In Figure 2A, the author claimed that “there was not transcriptome differences between with and without spike in”. This comparison is not rigorous. Importantly, the authors need to perform differential gene expression analysis (shown by Volcano plot) of (a) sorted 100% infected CD4 with HIV-seq versus (b) sorted 100% infected CD4 with regular 5’ CITEseq to see whether HIV spike in primers biased the cellular transcriptome. This is because HIV spike in primers may nonspecifically capture human transcriptome, or may eat into the regular cellular transcriptome capture. Starting with the viremic samples proves that the assay “sort of” works but without comparison to the ground truth.\n3. Figure 2C and 2D: without an uninfected donor control to set the threshold, the HIV+ cells identified in this plot cannot be distinguished from false positive detection. Again, a negative control and a threshold is needed. For example, “HIV transcript levels ranging from 1 to 1063 HIV reads per cell” – the field has identified that using one HIV RNA read to call HIV RNA+ cell yield false positive detection. Again, uninfected donor HIV-seq and regular CITEseq is needed – otherwise the data of high levels of HIV detection cannot be trusted.\n4. Please do not use “adjusted p value ≤ 0.05” – it has to be <0.05, not “less and equal”. There are a lot of “p <= 0.05” in the paper – please correct them all. This is fundamental statistical mistake.\n5. Figure 3E. How do you explain that pathway analysis of DEG in Figure 3E did not show cytolytic-related signatures? HIV RNA+ were also enriched in cluster 1 (~30%), what are the signatures for cluster 1?\n6. Figure 3F and 3G, Figure 4E, Figure S5. Please do not leave clusters as 1,2,3,4,5,6 without annotation. Please refer to major single cell papers in the immunology, cancer, or HIV field, such as Figure 2A-2B of PMID 31359002 or Figure 3A-3B of PMID 36536105. All clusters need to have dot plots or signature gene expression and heatmaps (of differentially expressed genes) to justify their cluster name.\n7. In Figure 3B, HIV RNA+ cells are enriched in the right part of the cluster which indicates cluster 1 from Figure 3F is a heterogenous cluster containing at least 2 distinct populations. The resolution for UMAP needs to be improved for further analysis and the inaccurate clustering results in biased conclusion. The same concern is also addressed on Figure 4E and 4F. The authors need to increase the resolution until biologically relevant CD4+ T cell clusters (supported by their cluster annotation by RNA and protein expression, shown by dot plots) are present and reflect human physiology. There are more than 6 different types of CD4s in the peripheral blood, but there are only 6 clusters without biological meanings.\n8. Figure 3C, 4C, and 5D showed that HIV+ cells were enriched in Tem. Can you visualize the Tem/Tcm/Ttm/Tn/CCR7+, CD37- populations on UMAP? Is the Tem population overlap with HIV RNA+ spots?\n9. Figure 4F. To claim HIV enrichment in IL7R and BCL2-high cells, using feature plot in UMAP is not rigorous. The authors should use violin plots of IL7R and BCL2 expression in HIV+ versus HIV- cells and perform statistical test (eg. Wilcox) with Benjamini-Hochberg corrected P values to rigorously justify this statement.\n10. Batch effect correction: Seurat default batch effect correction (integration) is known to have problems. For example, Figure S6 shows batch effects not corrected. Please use a different batch effect correction method (such as Harmony) and replot Figure S6 to see if there are still batch effects.\n11. Removing HIV reads from clustering: it’s great that the authors are aware that “T cell Receptor Alpha Variable (TRAV) and T cell Receptor Beta Variable (TRBV) genes were removed from the variable features used for clustering as they were driving the clustering in a donor-specific manner”. In fact, HIV leads to the same bias. Please remove HIV from the clustering process and re-define the clusters.\n12. Figure 1 – except for spiking in HIV, Figure 1B is generic from 10x platform that is already commercially available and should not be presented as Figure 1 as if it’s a novel design by the authors. Please remove Figure 1B.\nMinor comments:\n1. Table 1 is not in publication quality. Please change font size, remove grid lines (only horizontal lines)\n2. Please describe HIV mapping cut-off (based on uninfected donor results), package, version.\n3. Supplement table of differential gene expression – please remove it from Word file and use an Excel file instead so that the following figures are not blocked from reading.\n4. Supplement tables: gene names have to be italicized. Figures and text: HIV gene names (gag) should be all lower case and italicized. Please refer to genomics nomenclature guide in Wikipedia or genomics textbooks.\n5. Adjusted p-value should be used to analyze DEG and pathway analysis as assessing multiple comparisons on gene expression will increase the chance of false positive cases when using raw p-value.\n6. Can the author provide information/data to demonstrate the titration of spike-in HIV probes? What are the differences between spiking in 0.33 µM vs. 0.67 µM primers? This will be very helpful for the field.", "reviewer_2": "Frouard et al. describes an improved method called HIV-seq to identify HIV RNA+ cells using single-cell RNA-seq (scRNA-seq) in samples from people living with HIV. This approach builds upon previous literature to identify HIV RNA+ cells using the 10x Genomics platform. In this manuscript, the authors introduce to the default 10x single cell workflow a capture-based method to target HIV transcripts and enrich cells that are actively transcribing HIV using probes that are specifically targeting B subtype HIV. Currently, detection of HIV RNA+ cells is a major challenge in single-cell studies, particularly at ART when viral transcription is low, thus the authors’ methodology could be a significant contribution to aid this area of research. The expectation leading into the results section of the manuscript was that this was largely a methods paper, which follows logically from the main text of the very well-written introduction. However, despite this premise, the results moved very quickly from implementation and validation of the actual method and into a more superficial analysis of differential gene expression with overinterpretation of findings with biology of HIV persistence. Overall, the authors have not provided sufficient evidence to support their claim that this enrichment method increased detection of HIV RNA+ cells over the previous standard 10x method.\nMajor comments:\nGiven the small increase in frequency of HIV RNA+ cells, the lack of the total number of cells in the text and figure legends makes accurate interpretation of the HIV RNA+ cells challenging. The absolute numbers of HIV RNA+ and total CD4 T cells per participant with and without HIV-seq will determine if there is actually a substantial increase in cells with HIV and should be included to assess if there is actually a consequential increase in HIV RNA+ cells. This incremental change does not warrant the manuscript title.\nOn the flip side the authors use absolute numbers to compare HIV RNA without providing the total number of CD4 T cells assessed for each participant. Thus, statements such as in the discussion section page 16 lines 344-352 using absolute numbers are not comparable. For instance, if many, many CD4 T cells from a donor are sequenced, there is higher probability of finding HIV RNA+ cells than when fewer cells are characterized. This impacts even figures like 2D which show comparisons of number of HIV transcripts rather than proportions (are these differences significant). Figures 2A and 3A are missing the total (CD3+CD4+ T cells), as well as number of cells with HIV. Page 23 also does not provide the input cells per 10x well of a chip, and so it is difficult to evaluate the impact of HIV-seq from these experiments.\nIn the methods section CD4+ T cells were purified using negative selection and then PBMCs from the same donor was spiked back at a ratio of 5:100 in order to have non-CD4+ T cells to establish gates for defining HIV RNA+ cells. Clustering methods for scRNA-seq such as WNN are well-enough established that such lab approaches are unnecessary. Could this have led to inadvertent contamination with other cell types? There are also a number of Ig genes identified as DEGs in table S3 and S5. Why are immunoglobulin genes expressed in CD4 T cells? Were these variable genes used in the clustering step or is there cell contamination from non-CD4 T cells?\nThe authors rely heavily on exploratory analyses when they fail to detect significant DEGs between various conditions which emphasizes the problem of small sample sizes. Use of ‘nominal p-values isn’t necessarily a fatal issue, but here they are repeatedly used, and then downstream analyses are further applied onto these non-significant findings. This overall makes some of the results feel tenuous and less convincing.\nThe authors generated surface protein data (CITE-seq) – on only a very limited panel of antibodies (15). It’s unclear how or where that data is being used. The authors have also performed DEG analyses, but similar DEP or cluster analyses is missing or limited and thus CITE-seq analyses is biased and incomplete.\nHow was the cutoff for HIV low and high determined? Is this arbitrary, perhaps why a difference is not observed when stratifying by this category (Page 8)? This also drives home the fact we need more HIV RNA+ cells, not just increased reads from the same cells as shown in Fig 2C-D.\nThere seems to be a marked bias shown in Fig 2D indicating a 3’ bias. Can the authors improve this by including probes for the 3’ LTR region of the virus? Also given that there are many subtypes of HIV, can the authors expand this method to other non-B subtypes to make this approach more generalizable to other global subtypes?\nThe authors have validated previous findings shown during viremia that TEMs have the most number of HIV RNA+ cells and express cytolytic markers in a pool of >1000 HIV RNA+ cells. However, their key findings at the ART suppressed timepoint is much weaker as only 26 cells were detected with HIV RNA. Additional samples to increase cell numbers at the ART timepoint like the viremic will be essential to have meaningful results for comparisons and title is thus overstated.\nMinor points\nFurther, as seen in the supplementary tables, a lot of host transcripts are long-noncoding RNA and thus not polyadenylated are sequenced when using the default 10x platform and so authors should tone down their suggestions that only poly-A sequences are targeted.\nPage 7 lines 129-130: figure 2A simply shows cluster distributions on UMAPs. This evidence is insufficient to support the assertion of “no difference in global gene expression.”\nPages 8,9 lines 174-182: can you show where these subclusters fall on your UMAP? Would go a long way to helping contextualize the data shown in subfigure 3B.\nFigure 3F: pie charts are not ideal way to show data. It is very difficult to interpret changes in relative representation for clusters other than 1 and 2.\nAre there sufficient cells to do TCR analysis?\nWere HIV transcripts included while clustering?", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Point-by-point Response to Reviewers\nWe thank the reviewers for their thoughtful comments. Please find below our point-by-point\nresponse to the reviewers’ comments in blue italics, followed by our responses.\nReviewer #1:\nSummary and Major Issues\nIn this paper, Frouard et al. developed a new single-cell sequencing method, HIV-seq, to perform\ntargeted enrichment of HIV transcripts and improve the efficiency for detecting the rare HIV RNA+\ncells from PWH. The authors incorporated this method with single-cell CITE-seq and profiled\nclinical samples from PWH before and after ART treatment. By comparing the phenotypes during\nviremia and ART, the author described the signatures of HIV-infected cells during viremia and\nafter viral suppression. The major strength of this study is the adding spike-in primers before cells\nare loaded into the 10x Chromium machine for targeted enrichment increased the capture of HIV\nRNA+ cells which increased HIV capture rate. The iterative adjustment of spike-in primers is a\ngood endeavor optimizing this assay.\nWe thank the reviewer for their thoughtful and constructive feedback. We appreciate their\nrecognition of the strengths of our study, particularly our spike-in of HIV primers to enhance the\ncapture of HIV RNA+ cells.\nThe major weakness is the lack of demonstration of sensitivity, specificity, and how these spike-\nin primers biased cellular transcriptome.\nAs detailed below, we have conducted additional analyses to further demonstrate the sensitivity\nand specificity of our approach, and to confirm that the spike-in primers do not bias the cellular\ntranscriptome.\nThe lack of threshold of setting HIV RNA+ reads make the data not trustable, and thus an\nuninfected donor control is needed to set the appropriate threshold and to avoid false positive\ncalling of HIV+ cells.\nAs detailed below, we have addressed this concern by performing new experiments using\nuninfected donor controls. These experiments allowed us to establish an appropriate threshold\nand minimize the risk of false-positive classification of HIV+ cells.\nThe cell cluster annotation and differential gene expression analysis are not rigorous and should\nbe improved.\nAs detailed below, we have now provided additional information and clarification on the clustering\nand differential gene expression analysis approaches we had used, which had implemented\nrigorous statistical methodologies.\nHaving p <= 0.05, not p <0.05, is a fundamental statistical flaw.\nWe have now clarified that all our “p <= 0.05” statements were in fact p < 0.05, and we changed\nthe manuscript accordingly.\nIf well addressed, this study can bring innovation and value for the field.\nMajor comments:\n1. For any new single-cell assay, a sensitivity and specificity assay, and how these spike-in\nprimers biased cellular transcriptome must be established first. Please primary CD4+ T cells\ninfected with HIV and sorted/enriched (to obtain 100% HIV infection) versus primary CD4+ T cells\nwith no infection, to perform HIV-seq versus regular 10x 5’ CITEseq [4 samples: (a) sorted 100% infected CD4 with HIV-seq, (b) sorted 100% infected CD4 with regular 5’ CITEseq, (c) uninfected CD4 with HIV-seq, (d) uninfected CD4 with regular 5’ CITEseq]. This benchmarking experiment\nis needed to address: (a) How should the authors set the threshold of HIV detection in uninfected\ncells, such as using HIV read count or UMI count as a threshold? (b) what’s the sensitivity\ndetecting HIV RNA+ cells increase by using HIV-seq?\nWe agree that establishing the sensitivity and specificity of any new single-cell assay is essential,\nas well as assessing how spike-in primers may bias the cellular transcriptome. However, we\ndisagree that the specific experimental conditions requested by the reviewer are necessary to\nachieve these goals. In particular, the suggested experiment would not work for multiple reasons.\nFirst, it is not possible to obtain a population of 100% infected primary CD4+ T cells, even using\na VSV-G pseudo typed virus, which circumvents the requirement for CD4 and co-receptor\n(CCR5/CXCR4) engagement1. Even if one were to use an X4-tropic reporter HIV virus and to use\nex vivo stimulated CD4+ T cells to increase cellular permissivity, at any given time only a few\npercent of the cells would be infected. Based on extensive experience in working with such\ninfected primary cells, bead-based enrichment would only marginally increase the proportions of\ninfected cells, while sorting of HIV-infected primary cells would be harsh on the infected cells\n(which are more fragile than uninfected cells) and result in cells of insufficient viability and quantity\nfor follow-up assays. Instead, the only way to obtain a pure population of 100% infected cells\nwould be to use a cell line model – e.g., one of the JLAT clones, as compared to uninfected Jurkat\ncells. However, all JLAT cells – even amongst a single clone – exhibit a range of HIV RNA\nexpression levels, as we have previously demonstrated at the single cell level2. Even sorting of\nGFP+ JLAT cells will result in a population of cells that will re-distribute into GFP+ and GFP- cells.\nThis would pose a problem in a benchmarking HIV-seq analysis, since one cannot tell whether\nlack of detection of HIV RNA in that cell is due to a sensitivity issue, or simply because that cell\nis not expressing HIV RNA. Furthermore, JLATs are very different from primary CD4+ T cells from\npatients, which is what HIV-seq is designed for.\nWe fully agree, however, that it is important to establish a clear threshold of HIV detection in\nuninfected cells. Hence, as detailed below in our response to comment #3, we have addressed\nthis question by performing HIV-seq on both cells from PWH as well as cells from people without\nHIV (PWOH), and we demonstrated that no HIV+ cells were detected in the cells from the PWOH.\nThese new data are presented in the revised manuscript (Lines 137-141 and in Fig.S3).\nTo address the increase in sensitivity of detection HIV RNA+ cells by using HIV-seq, we had\ndirectly compared CD4+ T cells from PWH processed using HIV-seq versus regular 5’ CITE-seq\n(Lines 132-139 and in Fig. 2B and C in original manuscript). These results had showed that HIV-\nseq significantly improved the sensitivity of detecting of HIV+ cells, and in the revised manuscript\nare presented in Lines 142-150 and Fig. 2B and C.\nThe issue of establishing how the spike-in primers may affect the cellular transcriptome is\naddressed in point #2 immediately below.\n2. Impact of spike in on transcriptome: In Figure 2A, the author claimed that “there was not\ntranscriptome differences between with and without spike in”. This comparison is not rigorous.\nImportantly, the authors need to perform differential gene expression analysis (shown by Volcano\nplot) of (a) sorted 100% infected CD4 with HIV-seq versus (b) sorted 100% infected CD4 with\nregular 5’ CITEseq to see whether HIV spike in primers biased the cellular transcriptome. This is\nbecause HIV spike in primers may nonspecifically capture human transcriptome, or may eat into\nthe regular cellular transcriptome capture. Starting with the viremic samples proves that the assay\n“sort of” works but without comparison to the ground truth.\nAs detailed above in Response #1, it is not possible to obtain a 100% pure population of cells\nwhere every cell is known to harbor HIV RNA. However, such a population is not needed to\naddress this point the reviewer makes, that it is very important to formally demonstrate that HIV-\nseq does not bias the cellular transcriptome. Although the UMAPs in Fig. 2A revealed no obvious\nchange in cellular transcriptome upon addition of the HIV capture sequences, we agree with the\nreviewer that a more quantitative assessment is warranted to determine whether any\ntranscriptional changes are elicited by HIV-seq. Hence, we performed differential gene expression\n(DEG) analysis on the two patient specimens that were analyzed in the absence vs. presence of\nHIV-seq. These new results are now presented as volcano plots (Fig. S2 in the revised\nmanuscript) and in new supplementary tables (Tables S3 and S4), and we have included them\nhere.\nIn total, implementation of HIV-seq resulted in only 4 DEGs in PID1052 and 6 in PID8027. The\nonly gene common to both, AL138963.4, is a long non-coding RNA3 with limited characterization.\nGiven that such low numbers of DEGs are in line with what would be expected with technical\nreplicates of scRNA-seq analysis of the same sample, we believe our original conclusion that\nthere is no significant global gene expression difference between the two experimental conditions\nholds. These additional data are now presented in Lines 129-133 of the revised manuscript.\n3. Figure 2C and 2D: without an uninfected donor control to set the threshold, the HIV+ cells\nidentified in this plot cannot be distinguished from false positive detection. Again, a negative\ncontrol and a threshold is needed. For example, “HIV transcript levels ranging from 1 to 1063 HIV\nreads per cell” – the field has identified that using one HIV RNA read to call HIV RNA+ cell yield\nfalse positive detection. Again, uninfected donor HIV-seq and regular CITEseq is needed –\notherwise the data of high levels of HIV detection cannot be trusted.\nWe thank the reviewer for this important suggestion, and we agree that having uninfected donors\nas a control to set the threshold of HIV+ cells is a better control than our original negative control\nthat consisted of comparing reads to those observed in CD8+ T cells, which are not permissive\nto HIV. Hence, we have now run side-by-side HIV-seq vs regular 5’ scRNA-seq on two uninfected\ndonors. These results are presented here:\nIt demonstrates that among all of the uninfected samples, not a single HIV transcript was\nidentified, either in the absence or presence of HIV-seq, in contrast to our infected donor PID8027,\nwhere HIV RNA+ cells were detected using both methods. Therefore, we are confident that our\nthreshold of 1 HIV read is appropriate for defining HIV+ cells. These new data, which have been\nadded to the manuscript as Fig.S3B, are now discussed in Lines 137-141 in the Results of the\nmanuscript.\n4. Please do not use “adjusted p value ≤ 0.05” – it has to be <0.05, not “less and equal”. There\nare a lot of “p <= 0.05” in the paper – please correct them all. This is fundamental statistical\nmistake.\nWe have clarified that all the p-values we had indicated were in fact p<0.05, not p≤0.05. This\ncorrection has been made in all figures and tables.\n5. Figure 3E. How do you explain that pathway analysis of DEG in Figure 3E did not show\ncytolytic-related signatures?\nWe would like to clarify that the conclusion about enrichment of HIV RNA+ cells among cytolytic\nCD4+ T cells was made in Fig. 3F and G, where we demonstrated that the majority of HIV RNA+\ncells belonged to cluster 2, which harbored a cytolytic signature. By contrast, Fig. 3E (and 3D)\ncorresponded to a different type of analysis, where all HIV RNA+ cells were analyzed as a single\npopulation of cells and compared to all HIV RNA- cells. Since HIV RNA+ cells are in clusters\nbesides cluster 2, and since cluster 2 includes both HIV RNA+ and HIV RNA- cells, it is not\nunexpected that the cytolytic signature did not come out of the total DEG analysis (Fig. 3D) nor\nthe pathway analysis of those DEGs (Fig 3E). In essence, the difference relates to the fact that\nFig. 3D and E consider HIV RNA+ cells as a uniform bulk population of cells, in contrast to Fig.\n3F and G, which leverage the heterogeneity of the HIV RNA+ cells.\nHIV RNA+ were also enriched in cluster 1 (~30%), what are the signatures for cluster 1?\nCluster 1 cells expressed high levels of SELL, TCF7, IL7R, RACK1, LEF1 and CCR7, suggesting\ntheir identity to be a central memory (Tcm) cells and naïve cells. The unique features of this\ncluster, as well as the other clusters analyzed, are now presented as a new supplementary figure,\nFig. S8. The persistence of HIV in this cluster is consistent with known persistence of HIV in the\nTcm compartment4, which is now discussed in the paper (Lines 268-270).\n6. Figure 3F and 3G, Figure 4E, Figure S5. Please do not leave clusters as 1,2,3,4,5,6 without\nannotation. Please refer to major single cell papers in the immunology, cancer, or HIV field, such\nas Figure 2A-2B of PMID 31359002 or Figure 3A-3B of PMID 36536105. All clusters need to have\ndot plots or signature gene expression and heatmaps (of differentially expressed genes) to justify\ntheir cluster name.\nWe have now annotated the six clusters as requested. Clusters were annotated based on DEGs\nexpressed in each of the clusters (using the FindClusters function in Seurat), as well as by\nmanually assessing for expression levels of markers of classic cellular subsets. In summary,\ncluster 1 was annotated as a mix of central memory and naïve CD4+ T cells (defined by high\nexpression of SELL, TCF7, IL7R, RACK1, LEF1, CCR7); clusters 2 and 5 as two types of\neffector/cytolytic CD4+ T cells (defined by their expression of GNLY, PRF1, GZMA, GZMB,\nGZMH, GZMM, LAMP1), with cluster 2 expressing higher levels of all cytotoxic markers as\ncompared to cluster 5; cluster 3 as proliferating CD4+ T cells (defined by high expression of\nMKI67, CDK1, PCNA, MCMC2, TYMS, RRM2, TOP2A, RACGAP1); cluster 4 as a small cluster\nof platelet-CD4+ T cell aggregates previously described in scRNA-seq datasets5 (and defined by\nhigh expression of GP9, PF4, PPBP); and cluster 6 as activated HLA-DR-expressing CD4 + T\ncells (defined by high expression of HLA-DRB1, HLA-DPA1, HLA-DPB1, CD74). We have now\nupdated Figures 3F and 4E and the manuscript text accordingly (Lines 228-234) with the\nannotations. We have also, as requested, displayed the signature gene expression profiles of\neach of the clusters in dot plot format, to justify the cluster names (new Fig. S8). The new Fig. S8\nis presented here:\n7. In Figure 3B, HIV RNA+ cells are enriched in the right part of the cluster which indicates cluster\n1 from Figure 3F is a heterogenous cluster containing at least 2 distinct populations. The\nresolution for UMAP needs to be improved for further analysis and the inaccurate clustering\nresults in biased conclusion.\nAlthough clustering can be performed in a subjective manner to increase the numbers of clusters\nusing user-defined criteria, the current standard in the scRNA-seq field is to use objective criteria\nto prevent over- or under-clustering6. Otherwise, with user-selected choice of the clustering\nresolution parameter, one could derive as many or as few clusters as desired without solid\nbiological basis. We opted against user-selected choice of clustering since it is currently advised\nagainst by statisticians in the scRNA-seq field, as it results in different types of clustering between\ndifferent scRNA-seq studies. Instead, we implemented an approach that is now common in\nscRNA-seq analysis (and detailed in our Methods section Lines 576-578 of the initial manuscript,\nLines 619-621 of the revised manuscript), where the criteria for choosing the correct resolution is\nbased objectively and statistically as one that resulted in reproducible clusters across biological\nsamples in the study. More specifically, by using a cross-validation approach, we evaluated the\nclustering of cells in a given biological sample based on predictions derived after holding out cells\nfrom this sample, and clustering the cells from all the remaining samples, as described in recent\nstudies7,8. We have now clarified in our revised manuscript that our clustering method uses a\ncross-validation approach to objectively establish the numbers of clusters in our dataset (Lines\n621-624).\nThe same concern is also addressed on Figure 4E and 4F. The authors need to increase the\nresolution until biologically relevant CD4+ T cell clusters (supported by their cluster annotation by\nRNA and protein expression, shown by dot plots) are present and reflect human physiology. There\nare more than 6 different types of CD4s in the peripheral blood, but there are only 6 clusters\nwithout biological meanings.\nWith regards to the numbers of types of CD4+ T cells in peripheral blood, or in fact the number of\ntypes of any subset of cells in any compartment, such a precise number cannot be established,\nas it depends on what specific pathway or readout one is looking at. Indeed, one could consider\nthere to be millions of types of CD4+ T cells in blood, as every individual CD4+ T cell has a\ndifferent transcriptional signature. We would also like to caution against making conclusions by\nlooking at UMAP and shapes of clusters in the UMAP. As discussed by Chari et al., generating\nUMAP figures for exploratory cell subset annotation can generate false clusters9. Although user-\ndetermined cluster resolution can be useful for analyses where one wants to make a point about\nspecific aspects of biology, we disagree with the suggestion to manually increase cluster\nresolution, and instead opt for the standard, objective clustering that is based on statistical rather\nthan subjective measures, as detailed above and now in Lines 619-624.\n8. Figure 3C, 4C, and 5D showed that HIV+ cells were enriched in Tem. Can you visualize the\nTem/Tcm/Ttm/Tn/CCR7+, CD37- populations on UMAP?\nWe have now added a visualization of the classical CD4+ T-cell subsets on the UMAP in the new\nFig. S7A, and this is presented here:\nIs the Tem population overlap with HIV RNA+ spots?\nBy comparing Fig. 3F with the new Fig. S7A, one can see that cluster 2, where HIV RNA+ cells\nare enriched, also exhibits the most overlap with Tem cells. This is assessed more quantitatively\nin the new Fig. S7B (and presented below), where we show using pie graphs that among all the\nclusters, cluster 2 is the most enriched for Tem cells. Hence, the Tem population does indeed\noverlap with the location of the HIV RNA+ cells.\n9. Figure 4F. To claim HIV enrichment in IL7R and BCL2-high cells, using feature plot in UMAP\nis not rigorous. The authors should use violin plots of IL7R and BCL2 expression in HIV+ versus\nHIV- cells and perform statistical test (eg. Wilcox) with Benjamini-Hochberg corrected P values to\nrigorously justify this statement.\nIn response to these comments, we first performed a statistical analysis of IL7R expression in\nHIV+ versus HIV- cells using Wilcoxon rank-sum tests with Benjamini-Hochberg correction:\nWe found that the data were more effectively visualized as bar graphs (top) rather than as violin\nplots (bottom), because of the relatively high numbers of HIV- cells forming the “tails”, making it\nvisually look like expression levels were higher among the HIV- cells when in fact it is in fact higher\namong the HIV+ cells. Hence, we have opted to display the comparison as bar graphs, which are\nnow presented as the new Fig. 4F. These results confirmed our initial report that IL7R expression\nis higher among HIV+ vs HIV- cells.\nWhen we performed a similar analysis for BCL2 expression, we found that the difference in BCL2\nexpression was not statistically different between HIV+ and HIV- cells:\nWe noted, however, that among all the HIV RNA+ cells, 44% expressed BCL2. Notably, all of\nthese BCL2+ HIV RNA+ cells expressed IL7R, and among the IL7R+ HIV RNA+ cells, 52.4%\nexpressed BCL2:\nOverall, these results that IL7R and BCL2 are commonly co-expressed, and that a substantial\nproportion of HIV RNA+ reservoir cells express both IL7R and BCL2, which was the original point\nwe had set out to make. These new analyses are now updated in the Figures 4F and G and\ndiscussed in Lines 270-275 of the revised manuscript.\n10. Batch effect correction: Seurat default batch effect correction (integration) is known to have\nproblems. For example, Figure S6 shows batch effects not corrected.\nWe would first like to point out that Seurat integration has been shown to be among the top three\nbatch correction methods in a widely cited benchmark study10. We had opted for Seurat\nintegration in our original analysis because it returns a batch-corrected count matrix that is then\nused for differential gene expression analysis. Harmony, another top-ranked batch correction\nmethod10, corrects low-dimensional representations (e.g., PCA, UMAP), but does not return a\nbatch-corrected count matrix, which complicates downstream analyses such as differential gene\nexpression analysis.\nPlease use a different batch effect correction method (such as Harmony) and replot Figure S6 to\nsee if there are still batch effects.\nIn response to the reviewer’s suggestion, we applied Harmony as an alternative batch correction\nmethod and re-plotted the original analyses performed in Figure S6. The analysis revealed that\nthe batch effect persisted both before (top row) after (bottom row) correction by Harmony:\nGiven this result and the fact that Harmony does not return a batch-corrected count matrix, we\nhave opted to maintain our original approach and have not included the Harmony-corrected data\nin the revised manuscript. However, we would be happy to include the plots above as a\nsupplementary figure if the reviewer feels this is important.\n11. Removing HIV reads from clustering: it’s great that the authors are aware that “T cell Receptor\nAlpha Variable (TRAV) and T cell Receptor Beta Variable (TRBV) genes were removed from the\nvariable features used for clustering as they were driving the clustering in a donor-specific\nmanner”. In fact, HIV leads to the same bias. Please remove HIV from the clustering process and\nre-define the clusters.\nWe would first like to clarify that our rationale for removing TCR genes from the clustering analysis\nwas to avoid donor-specific clustering, as TRAV and TRBV genes have donor-specific features.\nThis donor-specific clustering would not be expected for HIV RNA transcripts, which were\ndetected in all the donors (i.e., unlike TRAV/TRBV, they were not unique to a particular donor).\nAn issue separate from the donor-specific clusters is whether HIV transcripts drove the clustering\noverall (as opposed to donor-specific clusters). Multiple analyses suggested that they did not.\nFirst, using Receiver Operating Characteristic (ROC) analysis11, we found that HIV transcripts\nwere not among the top 20 genes driving the clustering. These data are now included as a table\nentitled “Top20genes_clusters” that can be found after the references listed at the end of this\ndocument. Second, we removed HIV transcripts and re-clustered the cells using the same\nprincipal components (PCs) and resolution. We then compared cluster membership before vs.\nafter HIV transcript removal. Visualization of the data by UMAPs demonstrated that overall cluster\ndistribution was similar with (left) vs. without (right) HIV transcript inclusion:\nThird, for a more quantitative analysis, we examined the confusion matrix comparing cluster\nmembership before vs. after HIV transcript removal. This approach enables assessment of the\nextent to which cells moved between clusters before vs. after HIV transcript removal.\nIn the confusion matrix plot above, the boxes show what proportion of cells remained in the same\ncluster before and after removing HIV transcripts. The brighter the color, the more cells stayed\nput. The observation that the yellow/green boxes are all along the diagonal (corresponding to the\ncondition where each cluster identity stayed put), while all other boxes are purple, demonstrated\nthat cluster membership overall did not overall change after HIV transcript removal. Finally, to\nformally quantitate cluster stability, we calculated the adjusted Rand index (ARI) score, a\nreflection of how similar two types of clustering are. An ARI score of 1 corresponds to 100%\nidentity. We found that the ARI score comparing clustering with vs. without HIV transcripts was\nvery close to 1 (0.96), demonstrating that removing HIV transcripts had minimal effect on\nclustering. Based on all these results, we have opted to keep the HIV transcripts in our clustering,\nas the most unbiased method of clustering and also because it did not markedly change cluster\ndistribution. We have not included in the revised manuscript the above analyses, but would be\nhappy to do so if the reviewer feels it is important.\n12. Figure 1 – except for spiking in HIV, Figure 1B is generic from 10x platform that is already\ncommercially available and should not be presented as Figure 1 as if it’s a novel design by the\nauthors. Please remove Figure 1B.\nIn response to this comment, we have now removed the portion of Figure 1 that referred to the\ncommercially available 10X platform pipeline while retaining the general workflow of HIV-seq.\nMinor comments:\n1. Table 1 is not in publication quality. Please change font size, remove grid lines (only horizontal\nlines)\nThe formatting issue has been addressed in the revised manuscript.\n2. Please describe HIV mapping cut-off (based on uninfected donor results), package, version.\nAs detailed above in the response to question #3, to establish the threshold for detecting HIV-\npositive cells and transcripts, we have now performed HIV-seq vs. the regular scRNA-seq pipeline\non samples from uninfected donors (new Fig S3B and Lines 137-141) to establish that that our\nthreshold of 1 HIV read is appropriate to define HIV+ cells, since no HIV transcripts were identified\nin any of the uninfected samples, irrespective of HIV-seq. Our mapping was done using\nCellRanger software version 6.0.2 (from 10X Genomics), as indicated in the original manuscript\nLine 565 (now Line 608 in revised manuscript).\n3. Supplement table of differential gene expression – please remove it from Word file and use an\nExcel file instead so that the following figures are not blocked from reading.\nThe supplementary tables of differential gene expression are now all presented in the excel file\nentitled “Supplementary_Tables”, as requested.\n4. Supplement tables: gene names have to be italicized. Figures and text: HIV gene names (gag)\nshould be all lower case and italicized. Please refer to genomics nomenclature guide in Wikipedia\nor genomics textbooks.\nWe have updated the supplementary tables, figures, and text to ensure gene names follow the\ncorrect formatting guidelines.\n5. Adjusted p-value should be used to analyze DEG and pathway analysis as assessing multiple\ncomparisons on gene expression will increase the chance of false positive cases when using raw\np-value.\nThroughout the manuscript, we primarily reported adjusted p-values – these p-values were\nimplemented in Fig S3D, Fig. 3D + Table S3, Fig. 4D + Table S4, and Fig. S7 + Table S6 of the\noriginal manuscript (now Fig S5D, Fig. 3D + Table S5, Fig. 4D + Table S6, and Fig. S10 + Table\nS8 in the revised manuscript). Unadjusted p-values were only used for two analyses that did not\nyield statistically significant DEGs after multiple testing correction - Fig. 5C + Table S5, and Fig.\n5E + Table S7 of the original manuscript (now Fig. 5C + Table S7, and Fig. 5E + Table S9 in the\nrevised manuscript), which were considered exploratory analyses. We had clearly indicated the\nexploratory nature of these two analyses by stating: “As an exploratory analysis, however, we\nassessed the DEGs with the lowest raw p-values to gain insights into potential cellular pathways”\nLines 265-267 (now Lines 288-291 in the revised manuscript) and “As an exploratory analysis,\nwe assessed the identities of the DEGs based on the raw p-values” Lines 299-300 (now Lines\n322-323 in the revised manuscript).\n6. Can the author provide information/data to demonstrate the titration of spike-in HIV probes?\nWhat are the differences between spiking in 0.33 µM vs. 0.67 µM primers? This will be very helpful\nfor the field.\nIn order to determine the optimal concentration for the spike-in HIV primers, we had tested\ndifferent concentrations across 3 experiments and on cells from both viremic and suppressed\nsamples. These raw data were not presented in the original manuscript, but given this request we\nnow include them in the new Fig. S1. In experiment 1, we had first compared the effect of adding\n0.67 µM primers versus no primers on the cDNA recovery yield. We observed a significant\nreduction in cDNA yield, suggesting that 0.67 µM was too high of a concentration. This finding\nwas true for both suppressed and viremic samples. Hence, we next reduced the amount of primer\nby half in two other experiments (to 0.33 µM) and found that, compared to no primers, this lower\ncondition did not impair cDNA yield. The cumulative results of all three experiments, which led us\nto choose the 0.33 µM concentration, are presented in the new Fig. S1 and are depicted below:\nReviewer #2\nFrouard et al. describes an improved method called HIV-seq to identify HIV RNA+ cells using\nsingle-cell RNA-seq (scRNA-seq) in samples from people living with HIV. This approach builds\nupon previous literature to identify HIV RNA+ cells using the 10x Genomics platform. In this\nmanuscript, the authors introduce to the default 10x single cell workflow a capture-based method\nto target HIV transcripts and enrich cells that are actively transcribing HIV using probes that are\nspecifically targeting B subtype HIV. Currently, detection of HIV RNA+ cells is a major challenge\nin single-cell studies, particularly at ART when viral transcription is low, thus the authors’\nmethodology could be a significant contribution to aid this area of research. The expectation\nleading into the results section of the manuscript was that this was largely a methods paper, which\nfollows logically from the main text of the very well-written introduction.\nWhile our study introduces HIV-seq as an improved method for detecting HIV RNA+ cells in\nscRNA-seq, we do not view it as a methods paper, as it provides the first extensive scRNA-seq\ncharacterization of HIV RNA+ cells from ART-suppressed people with HIV (PWH) and their direct\ncomparison to HIV RNA+ cells during viremia from the same individuals. Such comparisons could\nnot previously be performed due to insufficient numbers of HIV RNA+ cells that could be\nrecovered by scRNA-seq analysis of specimens from ART-suppressed PWH. We also note that\nthe characterization of 1,072 HIV RNA+ from viremic individuals is the largest such dataset to\ndate, which will be a useful resource for the field.\nHowever, despite this premise, the results moved very quickly from implementation and validation\nof the actual method and into a more superficial analysis of differential gene expression with\noverinterpretation of findings with biology of HIV persistence.\nWe are not clear exactly how the reviewer finds our differential gene expression analysis\n“superficial” nor which aspect of our findings he/she felt we were overinterpreting, but we would\nlike to point out that our gene expression analyses were all performed using rigorous statistical\nmethodologies. We used a stringent QC pipeline to ensure analysis of only high-quality single-\ncell data (e.g., removing cells with high mitochondrial gene expression (>15% mtDNA) and low\nfeature counts (<500)), and additionally removed doublets using the recently published approach\nDoubletFinder12. Random Forests and silhouette score-based assessment of clustering validity\nand subject-wise cross-validation approach were implemented to ensure a statistically-validated\nmethod of clustering. For gene expression analysis, we aggregated cell-level data into\npseudobulk samples grouped by donor before applying edgeR13,14 (a differential expression\ntesting method for RNA-seq datasets) which improves the power of differential expression\nanalysis by accounting for inter-individual variability, thereby reducing false positives13. Finally,\nwe had reported p-values that passed the 0.05 thresholds after adjustment for multiple correction,\nand log2 fold-change > 0.25, as used by others in the field15–17; in two instances only – where we\nclearly stated that we were performing an exploratory analysis – we reported results that were\nfound based on adjusted p-values, as detailed further below in the response to Comment #4. If\nthere are specific aspects of differential gene expression or specific instances of\noverinterpretation of findings that were not addressed in the comments that follow, we would be\nhappy to address them.\nOverall, the authors have not provided sufficient evidence to support their claim that this\nenrichment method increased detection of HIV RNA+ cells over the previous standard 10x\nmethod.\nAs detailed further in the responses to Comments #1 and #2 below, we had shown that HIV-seq\nincreases both the number of HIV RNA+ cells detected (by up to 44%) (Fig.2B/Lines 134-137 in\noriginal, Fig.2B/Lines 142-148 in revision) and the number of HIV transcripts per cell (with the\nmean reads per cell increasing from 22 to 44) (Fig.2C/Lines 137-139 in original, Fig.2C/Lines 148-\n150 in revision).\nMajor comments:\n1. Given the small increase in frequency of HIV RNA+ cells, the lack of the total number of cells\nin the text and figure legends makes accurate interpretation of the HIV RNA+ cells challenging.\nThe absolute numbers of HIV RNA+ and total CD4 + T cells per participant with and without HIV-\nseq will determine if there is actually a substantial increase in cells with HIV and should be\nincluded to assess if there is actually a consequential increase in HIV RNA+ cells. This\nincremental change does not warrant the manuscript title.\nWe had initially only reported the frequencies of HIV RNA+ cells with versus without HIV-seq, and\nin response to this request have now additionally added the absolute numbers of both HIV RNA+\ncells and total CD4 + T cells (Lines 142-148 in revised manuscript). We note that the manuscript\ntitle does not say anything about the magnitude of increase in detection of HIV RNA+ cells.\nInstead, it focuses on the notion that the host gene expression profiles are different among HIV\nRNA+ cells during viremia vs. during ART suppression. Our data strongly support this finding,\nirrespective of the magnitude of increase that HIV-seq enables in recovery of HIV RNA+ cells.\n2. On the flip side the authors use absolute numbers to compare HIV RNA without providing the\ntotal number of CD4 + T cells assessed for each participant.\nWe have now added the total number of CD4 + T cells assed in our study (Lines 166 and 256)\nand per participant (Figure legend Lines 994-996 and 1028-1030). Specifically, 85,667 CD4+ T\ncells from 4 viremic samples (13,035, 35,849, 17,723, and 19,060 CD4+ T cells for PID1052,\nPID8026, PID8027, and PID145, respectively) and 75,158 CD4+ T cells from 3 suppressed\nsamples (15,468, 42,618 and 17,072 CD4+ T cells for PID1052, PID8026 and PID8027,\nrespectively) were analyzed.\nThus, statements such as in the discussion section page 16 lines 344-352 using absolute\nnumbers are not comparable. For instance, if many, many CD4 + T cells from a donor are\nsequenced, there is higher probability of finding HIV RNA+ cells than when fewer cells are\ncharacterized.\nAs detailed above, our study had sequenced a total of 85,667 CD4+ T cells to identify 1,072 HIV\nRNA+ cells during viremia. We have now detailed in the aforementioned Discussion section how\nour having identified many more HIV RNA+ cells was not due to our having sequenced many\nmore cells than prior studies. Specifically, the Collora et al.15 study, which had identified 81 HIV\nRNA+ cells, had come from sequencing 52,473 CD4+ T cells. Comparing to our study, we had\nsequenced 1.6x more cells, but found 13.2x more HIV RNA+ cells, demonstrating that our finding\nmore HIV RNA+ cells was not because we had sequenced many more cells. The other two prior\nstudies that had used single-cell sequencing had identified 164 HIV RNA+ cells from sequencing\n95,872 PBMCs (in Geretz et al.17), and 256 HIV RNA+ cells from sequencing 25,778 memory\nCD4+ T cells (Wei et al.16). Because the exact numbers of CD4+ T cells these corresponded to\nwere not provided in these manuscripts, such a direct comparison as performed with the Collora\net al.15 study is not possible. Regardless, if we approximate that CD4+ T cells comprise up to 60%\nof PBMCs18, then we estimate that 57,523 CD4+ T cells were sequenced in the Geretz et al.17\nstudy, which equates to our having sequenced 1.5x more CD4+ T cells but having identified 6.5x\nmore HIV RNA+ cells. Likewise, for the Wei et al.16 study, if we assume that total CD4+ T cells in\nPBMCs are comprised of ~50% of memory CD4+ T cells19, then that study sequenced the\nequivalent of 51,556 total CD4+ T cells. This information equates to our having sequenced 1.6x\nmore CD4+ T cells but detected 4.2x more HIV RNA+ cells. Hence, these results demonstrate\nthat our recovery of HIV RNA+ cells was unlikely the result of sequencing more input cells. We\nhave now modified the manuscript to better reflect this point (Lines 375-385).\nThis impacts even figures like 2D which show comparisons of number of HIV transcripts rather\nthan proportions (are these differences significant).\nIn response to this comment, we have modified Fig. 2D to report HIV transcripts as a proportion\nof CD4 + T cells rather than as absolute numbers. As seen below, the original data reporting\nabsolute numbers of reads (left) look very similar to the updated figure showing normalized data\nreporting their proportions (right):\nFigures 2A and 3A are missing the total (CD3+CD4+ T cells), as well as number of cells with HIV.\nWe believe the reviewer was referring to Figures 3A and 4A since Figure 2A is a UMAP. We have\nnow added the total number of CD4+ T cells analyzed to these figure panels (Line 167 for Fig. 3A\nand Line 256 for Fig. 4A). For consistency, in other figure panels where in the original manuscript\nwe had reported the absolute numbers and/or frequencies of HIV RNA+ cells, we now also always\ninclude the absolute numbers of CD4+ T cells as well (Lines 994-996 and 1028-1030 in the\nrevised manuscript).\nPage 23 also does not provide the input cells per 10x well of a chip, and so it is difficult to evaluate\nthe impact of HIV-seq from these experiments.\nWe had detailed in the original manuscript that we had added 20 µl of cells at ~1,000 cells/µl into\nthe reaction (Lines 545-546 in original manuscript). This amount corresponds to ~20,000 input\ncells when one does the math. In the revised manuscript, we have now explicitly stated that we\nhad inputted ~20,000 cells (Lines 588).\n3. In the methods section CD4+ T cells were purified using negative selection and then PBMCs\nfrom the same donor was spiked back at a ratio of 5:100 in order to have non-CD4+ T cells to\nestablish gates for defining HIV RNA+ cells. Clustering methods for scRNA-seq such as WNN\nare well-enough established that such lab approaches are unnecessary. Could this have led to\ninadvertent contamination with other cell types?\nThe PBMC spike-in was included in order to provide a population of cells non-permissive to HIV\ninfection (e.g., CD8+ T cells) to serve as negative controls for HIV RNA+ cell identification. This\npopulation was only used in the initial analysis shown in the original Figure S1 (now Fig. S3).\nFor all downstream analyses, we focused on a rigorously gated population of CD4+ T cells, defined\nas CD19-/CD20-, CD45+, CD3+, and CD8- cells. Furthermore, all samples were subjected to\ndoublet cell cleanup step using DoubletFinder12, minimizing the likelihood of doublets between\nCD4+ T cells and other cell types. Given these stringent QC and upstream gating criteria,\ncontamination with other PBMC-derived cell types would be minimal.\nThere are also a number of Ig genes identified as DEGs in table S3 and S5. Why are\nimmunoglobulin genes expressed in CD4 + T cells?\nAlthough we had followed 10X Genomics guidelines for cell loading to keep the doublet rate low\n(by loading a maximum number of ~20,000 cells to keep the doublet rate below ~7.6%), and\nimplemented DoubletFinder as an additional method to bioinformatically exclude doublets, these\napproaches are not able to remove every last doublet. This result is expected and unavoidable in\nscRNA-seq experiments. We suspect that the presence of Ig genes in our DEG tables is due to\nsome residual doublets between CD4+ T cells and B cells. A second possibility is that sample\nprocessing led to some lysis of B cells, which released B cell-derived ambient RNA which got\nencapsulated with CD4+ T cells. Both scenarios are well-known technical limitations of single-cell\nRNA-seq platforms.\nWere these variable genes used in the clustering step or is there cell contamination from non-\nCD4 + T cells?\nSince we could not determine the exact source of Ig gene signals in our data, and because there\ncould be other non-Ig genes that resulted from similar kinds of contamination, to be unbiased as\npossible, we opted to leave all genes that passed QC in during the clustering steps, rather than\nselectively removing ones that did not “make sense”.\n4. The authors rely heavily on exploratory analyses when they fail to detect significant DEGs\nbetween various conditions which emphasizes the problem of small sample sizes. Use of ‘nominal\np-values isn’t necessarily a fatal issue, but here they are repeatedly used, and then downstream\nanalyses are further applied onto these non-significant findings. This overall makes some of the\nresults feel tenuous and less convincing.\nWe would like to point out that throughout the manuscript, we primarily reported adjusted p-values\n– these p-values were implemented in Fig S3D, Fig. 3D + Table S3, Fig. 4D + Table S4, and Fig.\nS7 + Table S6 of the original manuscript (now Fig S5D, Fig. 3D + Table S5, Fig. 4D + Table S6,\nand Fig. S10 + Table S8 in the revised manuscript). Unadjusted p-values were only used for two\nanalyses that did not yield statistically significant DEGs after multiple testing correction - Fig. 5C\n+ Table S5, and Fig. 5E + Table S7 of the original manuscript (now Fig. 5C + Table S7, and Fig.\n5E + Table S9 in the revised manuscript), which were considered exploratory analyses. We had\nclearly indicated the exploratory nature of these two analyses by stating: “As an exploratory\nanalysis, however, we assessed the DEGs with the lowest raw p-values to gain insights into\npotential cellular pathways” Lines 265-267 (now Lines 288-290 in the revised manuscript) and\n“As an exploratory analysis, we assessed the identities of the DEGs based on the raw p-values”\nLines 299-300 (now Lines 322-323 in the revised manuscript).\n5. The authors generated surface protein data (CITE-seq) – on only a very limited panel of\nantibodies (15). It’s unclear how or where that data is being used. The authors have also\nperformed DEG analyses, but similar DEP or cluster analyses is missing or limited and thus CITE-\nseq analyses is biased and incomplete.\nOur CITE-seq panel was intentionally designed with a limited number of antibodies, with the\nprimary goal to facilitate gating, rather than as an unbiased way of gene expression analysis\nbased on protein expression. In particular, we used CD45 to identify immune cells; CD19, CD20,\nCD14, and CD8 to exclude B cells, myeloid cells, and CD8+ T cells; CD3 and CD4 to identify\nCD4+ T cells; and CD45RA and CD45RO to distinguish between naïve and memory CD4+ T\ncells. This information was shown in Fig. S4 of the original manuscript (now Fig S6 of the revised\nmanuscript). In addition, we included CD197 (CCR7), CD27, CD62L, CD127, CD279 (PD1), and\nCD49d as additional phenotyping markers and markers to help identify some cellular subsets like\ncentral memory T cells, effector memory T cells, and regulatory T cells. We have now clarified\nthe purpose of these phenotyping markers in the text (Lines 564-568 in revised manuscript).\n6. How was the cutoff for HIV low and high determined? Is this arbitrary, perhaps why a difference\nis not observed when stratifying by this category (Page 8)? This also drives home the fact we\nneed more HIV RNA+ cells, not just increased reads from the same cells as shown in Fig 2C-D.\nThe cutoff for HIV low and high was based on the observed distribution of HIV RNA+ cells in our\ndata, since we identified what seemed like a bimodal distribution of two populations of high vs low\ncells in the 2D dot plots shown in Fig. S3A and B of the original manuscript. The threshold of 50\nreads was chosen and applied across all samples, as this seemed to be the threshold that\ndistinguished the HIV RNA+ cells with lower vs. higher levels of HIV transcripts. Because a\nhistogram depiction of the data may better illustrate the separation between the HIV vs. HIV\nlow high\npopulations, we have now added that visualization as a new plot in Fig. S5B. That visualization is\ncopied below, where it is apparent, particularly when looking at the total population of HIV RNA+\ncells, that there is a bimodal distribution of HIV RNA+ cells:\nWith regards to needing more HIV RNA+ cells and not just increased reads from the same cells,\nwe note that HIV-seq does in fact increase the numbers of recovered HIV RNA+ cells as\ncompared to standard 5’ 10X sequencing. This result was demonstrated in Figure 2B in our\noriginal manuscript. We also note, as already discussed above, that our study reports the highest\nnumber of HIV RNA+ cells analyzed to date.\n7. There seems to be a marked bias shown in Fig 2D indicating a 3’ bias. Can the authors improve\nthis by including probes for the 3’ LTR region of the virus? Also given that there are many subtypes\nof HIV, can the authors expand this method to other non-B subtypes to make this approach more\ngeneralizable to other global subtypes?\nWe agree that additional probes targeting the 3’ LTR region would be of interest to test for\nimprovement, as would expanding the probes to enable analysis of other non-B subtypes and\ntesting specimens obtained from non-B regions of the world. However, these aspects are well\nbeyond the scope of the current study.\n8. The authors have validated previous findings shown during viremia that TEMs have the most\nnumber of HIV RNA+ cells and express cytolytic markers in a pool of >1000 HIV RNA+ cells.\nHowever, their key findings at the ART suppressed timepoint is much weaker as only 26 cells\nwere detected with HIV RNA. Additional samples to increase cell numbers at the ART timepoint\nlike the viremic will be essential to have meaningful results for comparisons and title is thus\noverstated.\nWe would like to point out that we detected the highest number of HIV RNA+ cells reported to\ndate in the context of ART suppression, and to date no published studies have compared HIV\nRNA+ cells in the context of ART suppression and viremia. One is limited by the low frequencies\nof HIV RNA+ cells that are present during ART suppression, and even if one were to spend years\nadding a lot more donors to the analysis (which is well beyond the scope of this study), the\nnumbers of HIV RNA+ cells would still overall be low due to the nature of the HIV reservoir.\nMinor points:\n1. Further, as seen in the supplementary tables, a lot of host transcripts are long-noncoding RNA\nand thus not polyadenylated are sequenced when using the default 10x platform and so authors\nshould tone down their suggestions that only poly-A sequences are targeted.\nWhile the 10X platform targets polyadenylated transcripts, it is known that the system is not\ncompletely clean and can capture some non-polyadenylated transcripts, likely through\nnonspecific priming or self-priming. We have now added a sentence stating this observation\n(Lines 106-107).\n2. Page 7 lines 129-130: figure 2A simply shows cluster distributions on UMAPs. This evidence\nis insufficient to support the assertion of “no difference in global gene expression.”\nWe agree with the reviewer that a more quantitative assessment of whether any transcriptional\nchanges is elicited by HIV-seq is warranted. Hence, we performed differential gene expression\n(DEG) analysis on two patient specimens that were analyzed in the absence vs. presence of\nHIV-seq. These new results are now presented as volcano plots (Fig. S2 in the revised\nmanuscript) and in new supplementary tables (Tables S3 and S4), and are depicted below:\nIn total, implementation of HIV-seq resulted in only 4 DEGs in PID1052 and 6 in PID8027. The\nonly gene common to both, AL138963.4, is a long non-coding RNA3 with limited characterization.\nGiven that such low numbers of DEGs are in line with what would be expected with technical\nreplicates of scRNA-seq analysis of the same sample, we believe our original conclusion of there\nbeing no significant global gene expression difference between the two experimental conditions\nholds. These additional data are now discussed in the revised manuscript (Lines 129-133).\n3. Pages 8,9 lines 174-182: can you show where these subclusters fall on your UMAP? Would go\na long way to helping contextualize the data shown in subfigure 3B.\nWe have now added a visualization of the classical CD4+ T-cell subsets on the UMAP in Fig.\nS7A.\nThis visualization demonstrates that the majority of HIV RNA+ cells overlap with the Tem\npopulation, supporting our statement that HIV+ cells are enriched in Tem cells.\n4. Figure 3F: pie charts are not ideal way to show data. It is very difficult to interpret changes in\nrelative representation for clusters other than 1 and 2.\nWe have now replaced the pie charts with stacked bar charts to indicate percentage distribution,\nsimilar to panel C, as a clearer representation of the data.\n5. Are there sufficient cells to do TCR analysis?\nUnfortunately, VDJ libraries were not generated for this study, so TCR analysis could not be\nperformed. As mentioned in the discussion of our initial manuscript (Line 468-471, Lines 501-504\nin revision), VDJ sequencing is compatible with the HIV-seq platform and should be integrated in\nfuture studies.\n6. Were HIV transcripts included while clustering?\nHIV transcripts were included in the clustering. The only genes removed were TCR genes, as\nthey caused donor-specific clustering driven by TCR expression. As detailed in Response #11 to\nReviewer #1, inclusion of HIV transcripts in the clustering did not alter cluster assignment nor\ndistribution.\nReviewer #3\n(Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nReferences:\n1. Bartz, S. R. & Vodicka, M. A. Production of high-titer human immunodeficiency virus type 1\npseudotyped with vesicular stomatitis virus glycoprotein. Methods 12, 337–342 (1997).\n2. Telwatte, S. et al. Heterogeneity in HIV and cellular transcription profiles in cell line models of\nlatent and productive infection: implications for HIV latency. Retrovirology 16, 32 (2019).\n3. Touil, H. et al. A structured evaluation of cryopreservation in generating single-cell\ntranscriptomes from cerebrospinal fluid. Cell Rep Methods 3, 100533 (2023).\n4. Chomont, N. et al. HIV reservoir size and persistence are driven by T cell survival and\nhomeostatic proliferation. Nat Med 15, 893–900 (2009).\n5. Dai, X.-P. et al. Increased Platelet-CD4+ T Cell Aggregates Are Correlated With HIV-1\nPermissiveness and CD4+ T Cell Loss. Front. Immunol. 12, (2021).\n6. Kiselev, V. Y., Andrews, T. S. & Hemberg, M. Challenges in unsupervised clustering of single-\ncell RNA-seq data. Nat Rev Genet 20, 273–282 (2019).\n7. George, A. F. et al. Deep Phenotypic Analysis of Blood and Lymphoid T and NK Cells From\nHIV+ Controllers and ART-Suppressed Individuals. Front Immunol 13, 803417 (2022).\n8. Yin, K. et al. Long COVID manifests with T cell dysregulation, inflammation and an\nuncoordinated adaptive immune response to SARS-CoV-2. Nat Immunol 25, 218–225 (2024).\n9. Chari, T. & Pachter, L. The specious art of single-cell genomics. PLOS Computational Biology\n19, e1011288 (2023).\n10. Tran, H. T. N. et al. A benchmark of batch-effect correction methods for single-cell RNA\nsequencing data. Genome Biology 21, 12 (2020).\n11. Jaskowiak, P. A., Costa, I. G. & Campello, R. J. G. B. The area under the ROC curve as a\nmeasure of clustering quality. Data Min Knowl Disc 36, 1219–1245 (2022).\n12. McGinnis, C. S., Murrow, L. M. & Gartner, Z. J. DoubletFinder: Doublet Detection in Single-\nCell RNA Sequencing Data Using Artificial Nearest Neighbors. Cell Systems 8, 329-337.e4\n(2019).\n13. Chen, Y., Chen, L., Lun, A. T. L., Baldoni, P. L. & Smyth, G. K. edgeR v4: powerful differential\nanalysis of sequencing data with expanded functionality and improved support for small\ncounts and larger datasets. Nucleic Acids Research 53, gkaf018 (2025).\n14. Robinson, M. D., McCarthy, D. J. & Smyth, G. K. edgeR: a Bioconductor package for\ndifferential expression analysis of digital gene expression data. Bioinformatics 26, 139–140\n(2010).\n15. Collora, J. A. et al. Single-cell multiomics reveals persistence of HIV-1 in expanded cytotoxic\nT cell clones. Immunity 55, 1013-1031.e7 (2022).\n16. Wei, Y. et al. Single-cell epigenetic, transcriptional, and protein profiling of latent and active\nHIV-1 reservoir revealed that IKZF3 promotes HIV-1 persistence. Immunity 56, 2584-2601.e7\n(2023).\n17. Geretz, A. et al. Single-cell transcriptomics identifies prothymosin α restriction of HIV-1 in vivo.\nSci Transl Med 15, eadg0873 (2023).\n18. Peripheral Blood | Whole blood | Handbook | Miltenyi Biotec | USA.\nhttps://www.miltenyibiotec.com/US-en/support/macs-handbook/human-cells-and-\norgans/human-cell-sources/blood-human.html.\n19. Su, L. F., Kidd, B. A., Han, A., Kotzin, J. J. & Davis, M. M. Virus-specific CD4+ memory\nphenotype T cells are abundant in unexposed adults. Immunity 38, 373–383 (2013).\nTop 20 genes clusters - Table showing the top 20 genes defining each cluster, identified using the\nFindAllMarkers function. These genes represent the most differentially expressed features\ndistinguishing each cluster in our analysis.\nmyAUC avg_diff power avg_log2FC pct.1 pct.2 cluster gene\n0.193 -2.3348246 0.614 -5.0708663 0.141 0.683 0 NKG7\n0.201 -1.7716477 0.598 -3.4886879 0.22 0.719 0 CCL5\n0.791 0.43818525 0.582 0.64997967 1 0.999 0 RPL32\n0.79 0.73823885 0.58 1.19528325 0.993 0.654 0 LTB\n0.788 0.45908093 0.576 0.67657052 1 0.999 0 RPS12\n0.777 0.34892301 0.554 0.51643183 1 1 0 RPL28\n0.225 -1.6189556 0.55 -4.7931754 0.037 0.571 0 GZMH\n0.773 0.39486338 0.546 0.59193811 1 1 0 RPL11\n0.771 0.41437924 0.542 0.61109316 1 1 0 RPLP1\n0.771 0.43842553 0.542 0.65612048 1 0.999 0 RPS8\n0.77 0.38614109 0.54 0.57321055 1 0.999 0 RPL30\n0.769 0.41228878 0.538 0.60948678 1 1 0 RPL13\n0.769 0.38016739 0.538 0.57853135 1 0.999 0 RPS25\n0.765 0.37794708 0.53 0.56470468 1 0.999 0 RPS28\n0.765 0.43650916 0.53 0.6634181 1 0.999 0 RPS13\n0.235 -2.2757943 0.53 -5.1860133 0.105 0.59 0 GNLY\n0.763 0.39254818 0.526 0.58523634 1 0.998 0 RPL39\n0.761 0.36442845 0.522 0.55275665 1 0.999 0 RPL12\n0.239 -1.2019387 0.522 -3.2953726 0.146 0.606 0 CST7\n0.76 0.37238775 0.52 0.55794255 1 0.999 0 RPS23\n0.961 2.04460403 0.922 3.3499365 0.992 0.333 1 NKG7\n0.942 1.67720456 0.884 2.78113042 0.992 0.401 1 CCL5\n0.931 1.38402699 0.862 2.82552553 0.968 0.276 1 CST7\n0.93 1.6042996 0.86 3.06857948 0.955 0.2 1 GZMH\n0.091 -1.7298964 0.818 -3.3615096 0.318 0.921 1 LTB\n0.9 1.45793818 0.8 3.17481355 0.889 0.123 1 FGFBP2\n0.888 2.050575 0.776 3.41806489 0.89 0.269 1 GNLY\n0.882 1.4440429 0.764 3.23313722 0.852 0.119 1 GZMB\n0.879 1.00010397 0.758 2.15672145 0.934 0.368 1 EFHD2\n0.873 1.1727343 0.746 2.61037603 0.879 0.179 1 PRF1\n0.87 1.06539276 0.74 2.22083248 0.918 0.255 1 GZMA\n0.862 0.50598516 0.724 0.74835661 1 1 1 HLA-B\n0.859 0.59219136 0.718 0.88894118 1 1 1 HLA-C\n0.854 0.52402468 0.708 0.76652955 1 1 1 B2M\n0.844 0.72108267 0.688 1.14212626 1 0.97 1 SH3BGRL3\n0.826 1.03353964 0.652 2.37453436 0.83 0.303 1 CTSW\n0.821 0.71981885 0.642 1.32381237 0.97 0.79 1 HCST\n0.821 0.59369248 0.642 0.92497755 1 0.992 1 PFN1\n0.817 0.4426398 0.634 0.67058973 1 1 1 HLA-A\n0.188 -0.9265131 0.624 -2.6190354 0.182 0.741 1 TCF7\n0.079 -0.9846992 0.842 -1.4755021 0.999 1 2 RPL10\n0.083 -1.006314 0.834 -1.4962126 1 1 2 EEF1A1\n0.903 0.9402016 0.806 1.36245764 1 1 2 MALAT1\n0.101 -0.8325547 0.798 -1.274441 0.999 1 2 RPS3\n0.102 -0.8777045 0.796 -1.3238554 0.999 1 2 RPL30\n0.108 -0.8416553 0.784 -1.2882958 0.997 0.999 2 RPS14\n0.109 -0.9046846 0.782 -1.4185212 0.997 1 2 RPS13\n0.109 -0.8919467 0.782 -1.3349694 0.999 1 2 RPS27\n0.111 -0.7586369 0.778 -1.1759852 0.996 0.999 2 RPL26\n0.113 -0.8297812 0.774 -1.2627057 0.998 1 2 RPS28\n0.117 -0.9256606 0.766 -1.3816689 0.998 1 2 RPS12\n0.118 -0.8496518 0.764 -1.3345021 0.997 0.999 2 RPS4X\n0.118 -0.716584 0.764 -1.1271185 0.997 1 2 RPS7\n0.124 -0.7352427 0.752 -1.1170065 0.999 1 2 RPL19\n0.126 -0.8466805 0.748 -1.2632566 1 1 2 RPLP1\n0.872 0.79121813 0.744 1.17196003 1 0.999 2 MT-CO2\n0.13 -0.7557711 0.74 -1.1603904 0.998 1 2 RPL18A\n0.131 -0.7614308 0.738 -1.1637982 0.998 1 2 RPS27A\n0.134 -0.7521028 0.732 -1.1477035 0.999 1 2 RPL11\n0.138 -0.7240252 0.724 -1.1379 0.994 0.999 2 RPL35A\n0.98 2.29985632 0.96 5.50034457 0.968 0.031 3 TRBV20-1\n1 2.82415375 1 6.94136211 1 0.017 4 TRBV7-2\n0.999 2.41445434 0.998 6.38179835 1 0.011 5 TRBV5-1\n1 2.60033664 1 6.71921488 1 0.011 6 TRBV7-9\n1 2.38726432 1 6.42766259 1 0.009 7 TRBV3-1\n1 2.49412821 1 6.70608967 1 0.004 8 TRBV10-3\n1 2.58392136 1 6.82737705 1 0.005 9 TRBV12-3\n1 2.4313496 1 6.55747636 1 0.005 10 TRBV4-2\n0.851 1.33492206 0.702 2.09668939 0.995 0.979 11 GAPDH\n0.849 0.87494011 0.698 1.33647809 1 0.994 11 PFN1\n0.849 0.94916228 0.698 1.40130609 1 0.999 11 ACTB\n0.848 0.82266214 0.696 1.37226808 0.99 0.934 11 CORO1A\n0.153 -1.0057309 0.694 -1.4652689 0.988 1 11 MALAT1\n0.838 1.13179836 0.676 2.21725103 0.899 0.529 11 COTL1\n0.826 0.76453967 0.652 1.21133529 0.999 0.987 11 CFL1\n0.176 -0.8627951 0.648 -1.5187685 0.519 0.977 11 TXNIP\n0.823 1.0467313 0.646 1.65300311 0.989 0.964 11 ACTG1\n0.822 0.76551645 0.644 1.40478496 0.958 0.795 11 ARPC1B\n0.181 -0.6645367 0.638 -1.0099831 0.986 0.999 11 RPL34\n0.819 0.70122999 0.638 1.18757261 0.989 0.94 11 MYL6\n0.817 1.56780111 0.634 3.1211505 0.839 0.472 11 HMGB2\n0.816 0.96359613 0.632 2.20255219 0.836 0.437 11 H2AFV\n0.199 -0.798039 0.602 -1.4934095 0.49 0.94 11 TSC22D3\n0.8 0.96999467 0.6 1.67887492 0.969 0.891 11 HMGB1\n0.2 -0.675725 0.6 -1.2437166 0.602 0.982 11 PNRC1\n0.798 0.79754462 0.596 1.69465526 0.869 0.58 11 COX8A\n0.797 0.6773138 0.594 1.41343343 0.901 0.611 11 PPP1CA\n0.796 0.77154122 0.592 1.83685837 0.825 0.451 11 COX5A\n0.833 1.98715108 0.666 5.58035609 0.677 0.025 12 NRGN\n0.795 1.74946172 0.59 5.29218694 0.602 0.021 12 PPBP\n0.776 1.10433395 0.552 4.22500165 0.56 0.013 12 CAVIN2\n0.751 1.19282288 0.502 4.43816744 0.509 0.01 12 TUBB1\n0.747 1.01356534 0.494 4.01214474 0.504 0.015 12 PF4\n0.726 0.67118 0.452 3.34235889 0.456 0.005 12 GP1BB\n0.718 1.22921748 0.436 3.6513003 0.516 0.111 12 HIST1H2AC\n0.717 0.60391533 0.434 3.14175127 0.439 0.006 12 SPARC\n0.713 0.70478873 0.426 3.36384472 0.433 0.01 12 PRKAR2B\n0.707 0.69500493 0.414 3.20714262 0.431 0.021 12 CLU\n0.707 0.65652881 0.414 1.1512443 0.94 0.873 12 TAGLN2\n0.982 2.58148531 0.964 6.74933735 0.965 0.007 13 TRBV7-3\n0.978 2.86109906 0.956 4.44292608 0.999 0.859 14 CD74\n0.964 2.71692735 0.928 5.07854232 0.975 0.244 14 HLA-DRB1\n0.964 3.20280084 0.928 6.51754877 0.95 0.102 14 HLA-DRA\n0.916 1.75660682 0.832 5.55305819 0.832 0.004 14 SPI1\n0.91 3.51891255 0.82 7.85394355 0.829 0.038 14 CST3\n0.903 1.75353154 0.806 3.99502665 0.864 0.22 14 TYMP\n0.901 3.27208046 0.802 7.77844856 0.805 0.016 14 LYZ\n0.899 2.80734641 0.798 6.85398105 0.808 0.037 14 IFI30\n0.899 1.46393418 0.798 2.75337003 0.923 0.531 14 COTL1\n0.899 2.46111843 0.798 4.65103957 0.912 0.334 14 HLA-DPA1\n0.894 1.67504181 0.788 3.3556875 0.9 0.461 14 PSAP\n0.894 1.90737662 0.788 4.89750714 0.813 0.088 14 AIF1\n0.892 1.25993186 0.784 1.92492655 0.997 0.997 14 FTH1\n0.888 1.54412487 0.776 4.27656703 0.813 0.113 14 GRN\n0.882 1.85444593 0.764 5.17524027 0.776 0.024 14 FCER1G\n0.877 1.04170704 0.754 1.95800423 0.958 0.765 14 GABARAP\n0.876 2.27183152 0.752 5.00069411 0.777 0.061 14 TYROBP\n0.875 1.59050928 0.75 3.61561266 0.853 0.295 14 CTSS\n0.873 1.55471984 0.746 4.57923166 0.77 0.067 14 LST1\n0.872 2.04287751 0.744 3.98678385 0.884 0.352 14 HLA-DPB1\nReviewer #1 (Remarks to the Author):\nIt is a pity that instead of addressing concerns and improving rigor, the authors did not\nrespond sufficiently to reviewer suggestions. This manuscript, despite presenting an\nexciting and innovative spike-in strategy, unfortunately did not demonstrate the rigor\ncompared with other manuscripts at Nature Communications.\nWe had applied rigorous methods to analyze the scRNA-seq data, and had performed\nnew experiments to demonstrate lack of detection of HIV RNA upon HIV-seq analysis of\ncells from people without HIV. As detailed further below, we have ensured that high\nrigor was applied in all aspects of the study, including adding a third uninfected donor\nand including details on all the numbers of cells sequenced from HIV- and HIV+\nsamples, which were similar (Lines 126-128, 139-142, 147-152, 170-172, 263-265,\n1029-1031). Although the number of clusters is not an issue of rigor, we have also\nincreased the resolution of the clusters, as suggested by the reviewer (lines 235-249,\n275-293; Fig. 2F, S9, S10, 3E). We are not aware of any additional issues with rigor.\nPoint 1. Sorting pure populations of HIV+ samples is technically feasible (PMID:\n30282021). The authors did not attempt to address this point.\nThere appears to have been a mis-understanding, as we were not making the case that\none cannot use a FACS machine to sort HIV+ cells (we agree it can be done, we\nourselves have extensively sorted HIV-infected primary cells, e.g., PMIDs 39798087,\n28746881). Rather, we were making the case that the “100% HIV infection” that the\nreviewer stated is not feasible to achieve with primary cells, since the initial infection\nfrequencies are generally low (<10%) and the sort purity maxes out at around 95%, so\nthere will be uninfected cells in the mix. As one cannot know the exact numbers of true\ninfected cells in these models, it will be difficult to establish the “upper end” of what\ncorresponds to detection of 100% of the infected cells. It was because of this that we\nhad discussed the notion of using clonally-derived cell line models (e.g., JLATs) to\nestablish what would be a 100% detection rate.\nBut more importantly, and which we apologize for not having clarified in our original\nresponse, there is a critical problem with trying to demonstrate that HIV-seq improves\ndetection of HIV RNA by performing a head-to-head comparison of conventional 10X\nvs. HIV-seq on either cell lines or even sorted HIV-infected primary cells. The problem is\nthat the types of HIV transcripts from cell line or primary cell models of HIV infection do\nnot recapitulate the types of HIV transcripts from ART-suppressed people with HIV\n(PWH) – in particular, these in vitro models do not replicate the excess of non-\npolyadenylated HIV transcripts that are observed in vivo. We have previously compared\nthe levels and patterns of different HIV transcripts in unstimulated cells from ART-\nsuppressed PWH to those measured in in vitro models of HIV persistence, including cell\nlines (7 different latently HIV-infected cell lines and the productively HIV-infected 8E5\ncell line1) as well as four different primary cell models of latent HIV infection2. The\nblocks to HIV transcription were different in each cell line or primary cell model, but\nnone of them showed the vast excess of non-polyadenylated HIV transcripts (block to\nHIV transcriptional completion) that we observed in cells from PWH. Specifically, in cells\nfrom ART-suppressed PWH, the ratio of 5’ elongated to polyadenylated (completed)\nHIV transcripts is about 10%, whereas in all cell lines and primary cell models that we\nhave studied, this ratio is at least 60% and usually close to 100% 1,2. In these cell lines\nand primary cell models, most HIV transcripts that have extended past the 5’ LTR have\na polyA tail, and therefore would be expected to be captured by the conventional 10x\napproach using a poly-dT. Hence, one would not necessarily expect any advantage to\nusing capture primers located between the 5’ LTR and the polyA tail, as used in HIV-\nseq. By contrast, HIV-infected cells in PWH – including in the context of ART\nsuppression as well as during active viremia3,4 – include high proportions of incomplete\ntranscripts, consistent with our ability to improve detection of HIV transcripts when\nspiking in the HIV-seq primers targeting non-polyadenylation signals. In summary, our\nHIV-seq was specifically designed to increase the capture of non-polyadenylated HIV\ntranscripts, which are primarily observed in vivo but not in vitro, and the best test of our\nmethod is to use actual samples from PWH, as performed in this manuscript. This is the\nmost rigorous approach and the gold standard, and the only way to really test the\napproach. We apologize for not having been clearer on this point in our original rebuttal,\nwhen we had simply stated that the in vitro models do not recapitulate in vivo infected\ncells – we were referring specifically to the types of HIV transcripts that are present,\nwhich we hope we have now clarified.\nPoint 3. While this reviewer appreciates the use of uninfected cells as a control as\nrequested, what's the number of cell sequenced, and how is the number comparable\nwith the HIV+ samples? It is showing a lack of rigor.\nWe apologize for leaving out these details. Because the numbers of cells sequenced is\nstandard and matched the numbers sequenced for the HIV+ samples, we had left out\nthose details. But we agree that for completeness sake and full transparency, we should\nshare all these details. For the two uninfected donors we had analyzed by HIV-seq, we\nhad sequenced 7,260 and 5,998 CD4+ T cells, which is within the range of cell numbers\nsequenced from the HIV+ study participants (median of 7,011 cells sequenced per\nperson; range 4,566-14,556). Because the numbers of cells sequenced for one of the\nuninfected donors was slightly on the lower end relative to the HIV+ samples, we have\nnow sequenced cells from an additional (third) uninfected donor, from whom we\nanalyzed 8,753 cells. In the three uninfected donors we analyzed by HIV-seq (where we\nsequenced 7,260, 5,998, and 8,753 cells), we did not identify a single HIV RNA read,\nwhereas in all of the samples from the HIV+ donors (where similar numbers of cells\nwere sequenced for each donor) we identified HIV RNA+ cells. These data clearly\ndemonstrate that HIV-seq does not result in false positive reads. These cell numbers\nare now detailed in lines 139-42 and the revised figure legend of Figure S4, while cell\nnumbers from the HIV+ participants are listed in lines 126-128, 147-152, 170-172, 263-\n265, and 1029-1031.\nPoint 7. For CD4 T cell resolution: the authors stated to stay with current resolution of 6\nclusters without minimal biological annotations. In other existing publications, there are\n10 clusters in PMID 35320704, 14 clusters in PMID 36070690, and 15 clusters in PMID\n37922905. While it is not the number of cluster that matters - it is whether the authors\nexamined and dissected the biology. See PMID: 33879890, Nature Med - 10 clusters,\nincluding naive, central memory, effector memory, Th1, Th2, Th17, Treg, proliferating, in\nthe peripheral blood. Leaving the annotation as cluster numbers without resolution to\nthese essential peripheral blood CD4 cell types does not bring insight to the field.\nPlease at least follow cell type annotations by PMID: 33879890.\nWe appreciate and thank the reviewer for bringing up other studies that had used\nbiological knowledge to justify increasing the numbers of clusters. Indeed, in the field\nthere is currently lack of consensus as to how best to define cluster resolution. While we\nstill think it is important to implement objective, statistics-based methods to define\ncluster resolution as suggested by others5, and as presented recently at the\nbioinformatics-based Systems for Molecular Biology (ISMB) conference6, and as we\nrecently implemented in other high-parameter single-cell sequencing studies7,8, we\nagree that bringing in biological knowledge is also important to establish biologically\nmeaningful clusters.\nTo bring both approaches together, we first implemented our statistics-based clustering\nmethod, and then selected the two largest clusters (clusters 1 and 2) that drove the\nunderclustering for further subclustering using the same optimized resolution applied in\nthe statistics-based approach. We then considered each of the subclusters as clusters,\nand added them to the remaining clusters that we deemed were not under-clustered\n(clusters 3-6 in original annotation). This approach resulted in a total of 13 clusters.\nWhen manually annotated, we could break down the clusters into the following subsets\n(note cluster numbers have been re-assigned, as there are now 13 clusters):\nCluster 1: Long-lived CD4+ T cells\nCluster 2: Tcm\nCluster 3: Non-cytotoxic Tem\nCluster 4: ISG+ long-lived CD4+ T cells\nCluster 5: GZMK+ Ttm\nCluster 6: Cytotoxic Tem\nCluster 7: KLRD1+ cytotoxic CD4+ T cells\nCluster 8: Tem\nCluster 9: Long-lived cytotoxic CD4+ T cells\nCluster 10: ISG+ Proliferating CD4+ T cells\nCluster 11: Platelet–CD4+ T cell aggregates\nCluster 12: Other CD4+ T cells\nCluster 13: Activated CD4+ T cells\nWe have incorporated the new clustering and cluster-distribution analyses into the\nResults (lines 235-249, 275-293), main figures (revised Fig. 2F and Fig. 3E), and\nsupplementary figures (Fig. S9, S10) as well as the new Supplementary Table 6, and\nwe revised the Methods and figure legends accordingly. Importantly, the main overall\nconclusions of our study remain unchanged: during viremia, HIV RNA+ cells\npredominantly reside in CD4+ T cells with a cytotoxic (high granzymes A/B/H/M,\nperforin, granulysin) and Th1 signature, whereas during ART suppression., they are\nfound not in cytotoxic CD4+ T cells but rather predominantly in a long-lived population of\nmemory CD4+ T cells. With the updated clustering, however, we were able to gain\nadditional insights into the features of HIV RNA+ cells during ART-suppression, in that\nwe found enrichment of HIV RNA+ cells in a second cluster that, in addition to having\nfeatures of stemness, displays high expression of MALAT1, a long non-coding RNA\npreviously shown to promote HIV transcription and infection9. These new findings have\nbeen added to the Results section of the revised manuscript (Lines 285-293).\nPoint 12: We requested removal of Figure 1B, which is a commercial 10x Genomics\nillustration - but Figure 1B remains in this version. The current Figure 1 remains\nextremely weak. I have never seen any publication with such a weak Figure 1 - probe\ndesign, that's it, no data, as a major Figure.\nIn our first revision, we removed the portion of Figure 1B that illustrated the\ncommercially available 10X platform pipeline while retaining the workflow unique to HIV-\nseq. We note that our revised Figure 1B from the previous submission was not a 10X\nGenomics illustration; it had showed an image of a 10X controller and an illustration\ncompletely created by us (not derived from any 10X Genomics literature) showing the\nrelationship of our HIV capture sequences to the PCR handle and to their binding\ntranscripts. In response to the reviewer’s request that this not be a main figure,\nhowever, we have now moved this schematic we created to a new supplementary figure\n(Fig. S2). In addition, we have merged the original Figure 1A with the previous Figure 2.\nThese panels now appear together as the new, consolidated Figure 1 in the revised\nmanuscript.\nReviewer #2 (Remarks to the Author):\nThe authors now provide the frequency of HIV RNA+ cells and the total number of CD4\nT cells assessed for each donor and experiment. From these numbers a z-test\nperformed by this reviewer, found no significant differences (P>0.05) for HIV RNA\nproportions with and without HIV-SEQ for all comparisons. It will be important to justify\nmore precisely than is currently provided, on what basis the inclusion of HIV-1 probes is\nclaimed to increase detection of HIV RNA+ cells.\nGiven the extremely high cost of scRNA-seq, the head-to-head comparison of HIV-seq\nvs. conventional 10x was performed on cells from a limited number of study\nparticipants. The z-test mentioned by the reviewer should not be applied to our data,\nsince as it assumes that the data follows a normal distribution and that the population\nvariance is known (not calculated from the data) or that the sample size is large (n>30) -\nconditions that are not met in our study. While the number of study participants does not\nsupport formal statistical comparison of HIV RNA+ cell frequencies, we found that the\nfrequency of HIV RNA+ cells increased by 26% and 44.7% in the two donors when\nusing HIV-seq (Lines 147-152), which supports the idea HIV-seq improves detection of\nHIV-infected cells. To further demonstrate that the HIV capture primers helped to\nincrease the detection of HIV RNA, we used the Mann-Whitney test to compare the\nnumber of HIV reads per cell detected with HIV-seq compared to the conventional 10x.\nWe found that HIV-seq doubled the mean number of reads per cell, and the difference\nwas very statistically significant (P<0.01; Fig. 1D in latest version of manuscript),\ndemonstrating with scientific rigor that HIV-seq enhances detection of HIV transcripts.\nWe now acknowledge that the small number of study participants did not support formal\nstatistical comparison of HIV RNA+ cell frequencies (lines 152-153), and we have\nrevised the text of the manuscript to emphasize that it is the number of reads per cell\nthat is statistically significant (lines 94, 152-154, 162-164, 389-390, 395).\nConsequently, although this paper is the second scRNA-seq study with >1000 HIV\nRNA+ cells detected, the authors have not currently demonstrated that this is because\nof their method.\nTo our knowledge, the only published studies specifically analyzing HIV RNA+ cells\nfrom people with HIV by scRNA-seq during viremia are the three we cited in the\nDiscussion (Collora et al. [25], Geretz et al. [32], and Wei et al. [51]). These studies\nreported a maximum of 256 HIV RNA+ cells (81 cells, 164 cells [with one highly viremic\ndonor further sequenced to gain an additional 223 HIV RNA+ cells] and 256 cells,\nrespectively).\nIt is possible that the >1,000 HIV RNA+ the reviewer referred to is from Figure 4 of the\nGeretz et al (reference 32) study, which was pointed out to the editor. In that figure,\nthere were numbers above the bar graphs corresponding to individual participants, one\nof which had a number greater than 1,000 under the category of “vRNA+”. While we\ncould not find definitions for what the numbers above the bars refer to within the Results\nor Figure Legends of the paper, given that this figure was brought up, we presume they\nrefer to numbers of HIV RNA+ cells that were identified. Going with this assumption, we\nhave now acknowledged in our revised Discussion that one prior study has detected\n>1,000 HIV RNA+ cells during viremia (Lines 401-403). We would like to point out that\nwe never made the claim that it is only because of HIV-seq that one can analyze >1,000\nHIV RNA+ cells. It is simply a method we present of interest for the field that increases\nthe efficiency of HIV RNA detection, and which we applied on cells from PWH and\nwhose datasets we then further analyzed. We have now added a sentence in lines 406-\n407 acknowledging that the high number of HIV RNA+ cells identified may reflect\ndifferences in the underlying frequency of HIV RNA+ cells as well as differences in HIV\nRNA capture.\nDetecting further HIV RNA+ cells at the viremic timepoint has only added incremental\ninsight into HIV RNA+ cells.\nWith due respect, and as detailed in our manuscript, we do not believe our analysis of\nthe viremic timepoint only added incremental insight into HIV RNA+ cells. New insights\nrelating to HIV RNA+ cells during viremia included:\n1) reduced expression of S100A8, S100A9, APOBEC3A, and SERPINA1 in HIV RNA+\ncells relative to HIV RNA- cells (lines 212-219, 437-452)\n2) upregulation of CXCR6 expression and an enhanced chemokine signaling pathway in\nHIV RNA+ cells compared to HIV RNA- cells (lines 224-225, 229-231)\n3) upregulation of both NFAT and PKC pathways in HIV RNA+ cells compared to HIV\nRNA- cells (lines 226-229)\nOur data also confirmed prior findings on HIV RNA+ during viremia (also important for\nthe field are confirmatory studies in different cohorts), in particular the demonstration\nthat:\n1) HIV RNA+ cells exhibit a heterogeneous phenotype but are enriched among Tem\ncells\n2) HIV RNA+ cells are enriched among cells with a cytotoxic/Th1 phenotype\nIn addition, we report new insights into the mapping of HIV RNA reads detected in\nsingle cells before and after ART (lines 160-162), the phenotypes and transcriptomes of\nHIV RNA+ cells during ART suppression (lines 260-301), and how they differ from HIV\nRNA+ cells during viremia (lines 341-385, 472-476, 477-526). Together with the new\nmethod used to increase detection of HIV RNA during scRNA-seq, our manuscript has\nadded much more than incremental knowledge. It should be noted that the description\nof our new method could comprise an entire manuscript by itself, yet we have combined\nit with a wealth of new data on the transcriptomes of cells in longitudinal samples from\nviremic and ART-suppressed people with HIV.\nNot using the ADT data in their differential analysis and generating TCR libraries are\nmissed opportunities.\nAs indicated in the prior rebuttal, we had used the ADT data primarily for cell\nphenotyping and gating rather than analysis of differential gene expression because the\nantibody panels were pre-selected and comprised a relatively small number of targets,\nas compared to the unbiased and full transcriptome data from scRNA-seq. However,\nwe appreciate the reviewer’s suggestion and agree that the limited protein expression\ndata can also be informative for differential gene expression analysis. Some markers—\nCD8, CD19, CD20, and CD14—were excluded from this analysis, as they were not\nrelevant for the CD4+ T cell population given that they were used to gate out CD8+ T\ncells, B cells, and monocytes. We performed differential protein expression analyses\nusing the remaining 11 markers and have now included these results in the revised\nmanuscript (Lines 200-207, 268-270, and 344-347; new Supplementary Figures S8,\nS11, and S14). We found that during viremia, HIV RNA+ cells express lower levels of\nCD45RA and higher levels of CD45RO compared to HIV RNA- cells across all donors,\nconsistent with our original reports that HIV RNA+ cells are depleted among naïve\nCD4+ T cells (which express more CD45RA) and enriched in memory CD4+ T cells\n(which express more CD45RO) (new Fig. S8). We also found some proteins that were\ndifferentially expressed in HIV RNA+ cells as compared to HIV RNA- cells, but only in a\nsubset of participants – these included CD279 (PD1; upregulated in HIV RNA+ cells in\n2/4 people), CD27 (downregulated in HIV RNA+ cells in 3/4 people), and CD49d\n(upregulated in 3/4 people; new Fig. S8). By contrast, during ART suppression we did\nnot find any ADT features consistently significantly differentially expressed among HIV\nRNA+ vs RNA- cells. However, comparing HIV RNA+ cells during viremia versus ART\nsuppression revealed lower CD4 expression in the context of active viremia, most likely\nreflecting either increased production or enhanced functionality of Nef protein, which is\nknown to downregulate cell-surface CD4 expression (new Fig. S14). These interesting\nnew data have now been added to our manuscript (Lines 344-347, new Supplementary\nFig. 14).\nRegarding TCR library generation, we agree this is a valuable tool, as it can inform\nabout clonal expansion dynamics and allow for prediction of antigen specificities. We\nhad acknowledged in the Discussion that multiplexing with TCR library generation is an\nimportant next step (Lines 528-532 revised manuscript), but note that these studies are\nwell beyond the current one, which was not focused on clonal expansion dynamics or\nantigen specificities.\nAlthough the authors sought to perform ‘a meaningful analysis of HIV RNA+ cells from\nART-suppressed PWH’ (line 99), they were only able to detect 25 HIV RNA+ cells at\nthis timepoint in 3 donors which did not allow for a rigorous analysis.\nAlthough 25 HIV RNA+ cells is not high, we note that this number is the highest\nanalyzed to date from blood of ART-suppressed people with HIV. Throughout the\nmanuscript, we had been careful not to overstate our findings. We explicitly indicated\nwhen analyses are exploratory and noted that such results should be interpreted with\ncaution. We also made sure that the results presented were rigorous and made use of\nappropriate statistical tests and data visualization approaches.\nFinally, and most importantly, reluctance to remove HIV while clustering,\nThe reviewer seems to be referring to a point brought up by the other reviewer\n(Reviewer #1, Point #11 from the prior revision). As detailed in our previous response\n(Response #11 to Reviewer #1), we had performed an analysis showing that removing\nHIV transcripts had minimal impact on clustering, given that 96% of the cells retained\ntheir original cluster identity. We had stated that this was our justification for not\nremoving HIV transcripts while clustering, and in our original response had offered to\ninclude these analyses if Reviewer #1 considered it important. Reviewer #1 made no\nfurther comments on the matter, so we assume Reviewer #1 does not feel it to be\nimportant.\nHowever, given Reviewer #2’s new concern on Reviewer #1’s prior comment, we have\nreplaced our initial clustering (performed with HIV genes included) with the version\ngenerated after excluding HIV genes from the clustering step (new Fig. 2F and Fig. 3E).\nThe Methods section has been updated accordingly (Lines 658) to reflect this change.\nAs expected, no conclusions are altered from this change, since 96% of the cells\nretained their original cluster identity.\nineffective batch correction,\nThis is another comment that was from the other reviewer (Reviewer #1, Point #10 from\nthe prior revision). Reviewer #1 had requested: “Please use a different batch effect\ncorrection method (such as Harmony) and replot Figure S6 to see if there are still batch\neffects.” We addressed this concern in our initial response (Response #10 to Reviewer\n#1) by showing that, even after applying the Harmony batch correction method as\nsuggested, donor-to-donor differences persisted. These results were presented as\nFigure S6 in our prior revision (copied again below for easy reference).\nThe differences observed with Harmony were similar to those seen with Seurat’s batch\ncorrection. As explained in our response to Reviewer #1, given that inter-donor variability\nis expected in primary human samples, we believe these differences reflect true biology\nrather than batch effects, and that no batch correction methods would be able to correct\nfor these biological differences. Reviewer #1 made no further comments on the matter.\nHowever, given this new concern from Reviewer #2 on the original comment from\nReviewer #1, we have replaced the plots in new Fig. S12 with the Harmony batch-\ncorrected version.\nand the absence of code to replicate their findings makes it challenging to evaluate the\nfindings and suggests a lack of rigor.\nWe note that the code for the quality control steps had already been shared previously,\nand the other steps involving code—such as clustering, differential expression analysis,\nand cell distribution-related statistics—were described in the Methods section, as is the\nstandard for scientific manuscripts.\nIn this second revision, we now show all code by providing a link to a Github repository\n(https://github.com/gladstone-institutes/HIV-seq-_Viremic_vs._ART) that has been\nmade available to the public.\nThis repository contains all scripts used in the manuscript:\n• Step0_QCsteps.v8.R – used for QC steps after CellRanger and for formatting the\ndata for SeqGeq software\n• Step1a.R and Step1b.R – used to format cell count data from SeqGeq into a\nSeurat object\n• Step2.R – used to analyze cell proportions for Figures 2C, 3C, 4B, and 4D\n• Steps3, 4, and 5 – used for normalization, visualization, and batch correction to\nprepare the data for visualization\n• Step6.R and Step7.R – used for cluster resolution optimization and cluster\ndefinition in Figures 2F and 3E\n• Step8.R – used for differential gene expression analysis in Figures 2D, 3D, 4C,\nand 4E\nThe Data Code and Availability Statement has been updated to include this link (lines\n727-733).\nOnly some data generated in the study was submitted to GEO.\nOur original GEO repository (GSE266455) only included the sequences generated from\nthe first version of the manuscript. We apologize that we had forgotten to upload the\nnew data generated for the previous round of revisions. To publicly share the new data\ngenerated for the revision, we have now created a new GEO repository (GSE305352).\nThis repository contains the two new uninfected donors analyzed as described in the\nprevious submission, in addition to one additional uninfected donor that we included\n(since one of the uninfected donors sequenced was on the lower end number of cells\nsequenced) as detailed in our response to Reviewer #1.\nThis sequence data can be accessed using the following links and tokens:\n----------\nTo review GEO accession GSE305352:\nGo to https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE305352\nEnter token ulwfmyesxxgtjgn into the box\n----------\nTo review GEO accession GSE266455:\nGo to https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE266455\nEnter token azodqesezbelnuj into the box\nThe Data Code and Availability Statement has been updated to include these links\n(lines 727-733).\nImproving existing methods is essential, but dependent on sufficient transparency for\nothers to validate results, which is not currently available here.\nWe agree that transparency is essential for validating scientific results. With the full set\nof codes now shared, the addition of the negative control data to the GEO repository,\nand our detailed responses to the reviewer’s remaining questions, we believe that we\nhave now provided sufficient transparency to enable others to fully assess and\nreproduce our findings.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports.\nThis is part of the Nature Communications initiative to facilitate training in peer review\nand to provide appropriate recognition for Early Career Researchers who co-review\nmanuscripts.\n---------------------------\n1. Telwatte, S. et al. Heterogeneity in HIV and cellular transcription profiles in cell line\nmodels of latent and productive infection: implications for HIV latency. Retrovirology 16, 32\n(2019).\n2. Moron-Lopez, S. et al. Human splice factors contribute to latent HIV infection in primary\ncell models and blood CD4+ T cells from ART-treated individuals. PLoS Pathog 16, e1009060\n(2020).\n3. Janssens, J. et al. Longitudinal changes in the transcriptionally active and intact HIV\nreservoir after starting ART during acute infection. J Virol 99, e0143124 (2025).\n4. Janssens, J. et al. Differential decreases in various HIV DNA regions and HIV transcripts\nafter ART initiation during chronic infection. J Virol 99, e0068325 (2025).\n5. Kiselev, V. Y., Andrews, T. S. & Hemberg, M. Challenges in unsupervised clustering of\nsingle-cell RNA-seq data. Nat Rev Genet 20, 273–282 (2019).\n6. Gill N & Thomas R. Optimizing Clustering Resolution for Multi-subject Single Cell\nStudies. Presented at Intelligent Systems for Molecular Biology (ISMB).\nhttps://github.com/gladstone-institutes/clustOpt (2025).\n7. George, A. F. et al. Deep Phenotypic Analysis of Blood and Lymphoid T and NK Cells\nFrom HIV+ Controllers and ART-Suppressed Individuals. Front Immunol 13, 803417 (2022).\n8. Yin, K. et al. Long COVID manifests with T cell dysregulation, inflammation and an\nuncoordinated adaptive immune response to SARS-CoV-2. Nat Immunol 25, 218–225 (2024).\n9. Wang, M.-R. et al. LncRNA MALAT1 Facilitates HIV-1 Replication by Upregulation of\nCHCHD2 and Downregulation of IFN-I Expression. Mol Cell Proteomics 24, 100997 (2025).\nPOINT-BY-POINT RESPONSE TO REVIEWERS\nReviewer #1 (Remarks to the Author):\nThis reviewer recognizes the effort that the authors attempted to revise the manuscript,\nalthough the quality of the bioinformatic analysis is not up to the current standard of the\nfield.\nTo our knowledge, every aspect of the bioinformatic analysis was performed according\nto the current gold standard of the field.\nThese are previously raised points that were not addressed.\nFigure 3A and 4E\nPlease annotate the clusters, instead of 1, 2, 3, 4, 5, 6.\nThis point was already addressed. We believe that this reviewer was looking at an\nolder version of the figures and not the most recent ones. The referenced number of\nclusters (6) and figure numbers (Fig. 3A and 4E) were from an older version of the\nmanuscript on which the reviewer requested cluster annotation. In the most\nrecent version of the manuscript, annotations were provided for clusters 1-13 in Fig. 2F\nand 3E, respectively .\nPlease integrate uninfected data into Figure 5.\nFigure 5A: one UMAP of batch-effect-corrected single cell data including data from\nviremic, suppressed, and uninfected donors.\nFigure 5B and 5D: Add uninfected donor data to the bar graph, as a reference to see\nhow HIV affects these T cells.\nFigure 5C: Add differential gene expression volcano plot of viremia vs uninfected,\nsuppressed vs uninfected\nFigure 5E: Add uninfected donor data to the bar graph, as a reference to see\nhow HIV affects these T cells.\nThis comment from Reviewer #1 appears to be a completely new request, and appears\nto be referencing an older version of the manuscript. In the latest version of the\nmanuscript, there is no Fig. 5 and main figures end at Fig. 4. While we could try to\nincorporate data from the uninfected donors, there is not enough space to include the\ndata in this figure, and it does not make sense scientifically. The 3 uninfected donors\nwere analyzed simply to demonstrate that the HIV-seq technique does not result in false\npositive detection of HIV transcripts. The figure referred to by the reviewer\ncompares paired samples from the same people before vs. after ART suppression.\nAdding in the 3 uninfected donors will result in a very mixed set of comparisons that\nwould be hard to interpret: we would be comparing longitudinal data from the same\ninfected people at two different timepoints (during viremia vs. ART suppression) but also\ntrying to compare data from those two timepoints cross-sectionally with data from 3\ndifferent HIV-uninfected donors. Perhaps even more importantly, the 3 uninfected\ndonors were in no way matched to the HIV+ individuals analyzed in this study (in age,\nsex, and other clinical aspects). Our study had not set out to compare HIV+ people to\nHIV- people, which would require a high \"n\" and a completely different study design.\nThis comparison should be reserved for a future study."}]

Decision letter