HIV silencing and cell survival signatures in infected T cell reservoirs – Nature

Study participants

Recruitment of study participants with HIV was performed in compliance with relevant ethical regulations under the IRB-approved SCOPE protocol (NCT00187512) at San Francisco General Hospital. Participants were enrolled from the SCOPE cohort on the basis of sample availability at the time of study, without use of sample size calculations, blinding or randomization. Demographic and clinical laboratory data were collected at San Francisco General Hospital and are reported in Supplementary Table 1. All of the participants provided informed consent before study. Prescreening of participant samples to ensure adequate numbers of HIV-DNA⁺ memory CD4 T cells for FIND-seq analysis was performed in parallel sample aliquots using fluorescence-assisted clonal amplification⁶³.

Cell lines

Jurkat human T cells (TIB-152, ATCC), HIV-DNA⁺ J-Lat full-length human T cells (clone 6.3, ARP-9846)⁶⁴ and Raji human B cells (CCL-86, ATCC) were cultured in Gibco RPMI Medium 1640 (Thermo Fisher Scientific, 11875093) with penicillin and streptomycin (Thermo Fisher Scientific, 15140122) and 10% fetal bovine serum (FBS). Mouse fibroblasts (NIH/3T3, CRL-1658, ATCC) were cultured in Dulbecco’s modified Eagle’s medium (DMEM) with penicillin and streptomycin (Thermo Fisher Scientific, 15140122) and 10% FBS. Before use, 3T3 cells were dissociated using 0.25% trypsin-EDTA (Thermo Fisher Scientific, 25200-072) and neutralized in DMEM with 10% FBS. Cell lines were used without authentication or mycoplasma contamination testing.

Fabrication of microfluidic devices

Standard photolithography techniques were used to fabricate microfluidic devices at the Harvard Medical School Microfluidics Facility. Silicon wafers were spin-coated with SU-8 2025/2050 photoresist (Kayaku Advanced Materials) and ultraviolet-patterned using a mask aligner. After developing, the wafers were baked overnight and used as master moulds for soft-lithography. In brief, the PDMS prepolymer and curing agent were mixed by hand at a ratio of 10 to 1 (Momentive, RTV615), degassed for at least 1 h, poured onto the mould and degassed until no bubbles remained. PDMS was baked overnight at 65 °C before holes were punched using a 0.75 mm biopsy punch and bonded to a glass slide (75 × 50 × 1.0 mm, Thermo Fisher Scientific, 12–550C) with a plasma bonder (Technics Plasma Etcher 500-II). Bonded devices were made hydrophobic with Aquapel with a 30 s contact time, flushed with HFE-7500, purged with air and baked for at least 1 h before use.

Cell line validation studies

Cells were washed twice with Hanks’ balanced salt solution (HBSS, no calcium, no magnesium, Thermo Fisher Scientific, 14170112) and then counted, mixed (mouse:human 1:1; J-Lat:Raji 1:100), and resuspended in HBSS containing 18% OptiPrep Density Gradient Medium (Sigma-Aldrich) for FIND-seq. For standard RNA-seq studies performed in parallel, aliquots of 5 × 10⁴ cells were lysed in RNAzol RT (Molecular Research Center) and stored at −80 °C until subsequent total RNA extraction according to the manufacturer’s instructions. Whole-transcriptome cDNA was then generated from total RNA by reverse transcription using 6 mM MgCl₂, 1 M betaine, 7.5% PEG-8000, 1 mM dNTP, 2 U µl⁻¹ Maxima H-minus reverse transcriptase (Thermo Fisher Scientific, EP0753), 0.5 U µl⁻¹ RNase inhibitor (Lucigen, NxGen) and 2 µM SMART TSO (AAGCAGTGGTATCAACGCAGAGTGAATrGrGrG). This cDNA was purified using AMPure XP beads (Beckman Coulter), and was then processed for WTA by PCR, with library preparation as previously described⁶⁵. FIND-seq sample processing and library preparation were performed as described below. The correlation between the DGE results from standard RNA-seq and FIND-seq was analysed using stat_cor (method = “pearson”) in R (v.4.1.0). The results from the J-Lat:Raji mixing study were compared with published transcriptomic signatures of CD4 T cells and B cells⁶⁶ using GSEA.

PBMC processing for FIND-seq

Approximately 20–30 million cryopreserved peripheral blood mononuclear cells (PBMCs) from each study participant were used for FIND-seq. Cryopreserved PBMC suspensions were thawed in a 37 °C water bath, washed in prewarmed RPMI with 10% FBS, and sedimented by centrifugation at 300 rpm (Sorvall Legend XT). Untouched memory CD4 T cells were then isolated by magnetic-column-based negative selection (Miltenyi, 130-091-893). Cells were counted manually with a haemocytometer using Trypan blue, and aliquots of 5 × 10⁴ cells were lysed and stored in RNAzol RT.

FIND-seq

FIND-seq was performed as described previously³⁰. In brief, four syringes were prepared for microfluidic cell encapsulation: lysis buffer, agarose, cells and oil. The lysis buffer consisted of 20 mM Tris-HCl pH 7.5, 1,000 mM LiCl, 1% LiDS, 10 mM EDTA, 10 mM DTT and 0.4 µg µl⁻¹ proteinase K. Conjugated agarose-dT was heated to 95 °C for 1 h before use and was kept heated throughout the run using a custom syringe heater. A 10 ml syringe was loaded with oil (Bio-Rad, 186–3005) for droplet generation. All of the syringes were connected to the microfluidic device using PE/2 tubing (Scientific Commodities, BB31695-PE/2). To make droplets, pumps were run at 600 μl h⁻¹ (cell mixture), 1,200 μl h⁻¹ (agarose), 600 μl h⁻¹ (lysis buffer), and 5,000 μl h⁻¹ (oil) using a bubble-triggered drop generator⁶⁷. Air was controlled to break the jet and generate 53–55 µm droplets. After lysis at 55 °C for 2 h, droplets were cooled at 4 °C overnight to allow agarose gelation. Solid agarose microspheres (beads) were removed from the oil using a drop-breaking procedure. All of the steps were performed at 4 °C to prevent dissociation of mRNA from the poly(T) oligonucleotides. The beads were removed from the oil and washed five times. For each wash, the beads were incubated in wash buffer for 5 min on ice, centrifuged at 4,700 rpm for 10 min and aspirated before the next wash. Beads were first washed in wash buffer 1 containing 20 mM Tris-HCl pH 7.5, 500 mM LiCl, 0.1% LiDS and 0.1 mM EDTA. Next, the beads were washed twice with wash buffer 2 containing 20 mM Tris-HCl pH 7.5 and 500 mM NaCl. Finally, the beads were washed twice in 5× reverse transcription buffer containing 250 mM Tris-HCl pH 8.3, 375 mM KCl, 15 mM MgCl₂ and 50 mM DTT and filtered with a 100 µm cell strainer. The beads were resuspended in reverse transcription master mix to a final concentration of 6 mM MgCl₂, 1 M betaine, 7.5% PEG-8000, 1 mM dNTP, 2 U µl⁻¹ Maxima H-minus reverse transcriptase (Thermo Fisher Scientific, EP0753), 0.5 U µl⁻¹ RNase inhibitor (Lucigen, NxGen) and 2 µM SMART TSO (AAGCAGTGGTATCAACGCAGAGTGAATrGrGrG). Reverse transcription was completed at 25 °C for 30 min, followed by 90 min at 42 °C. The tubes were mixed continuously with an inverter during all incubations. After reverse transcription, the beads were washed five times with 0.1% Pluronic in RNase/DNase-free water.

After reverse transcription, the cell occupancy of agarose beads was quantified by microscopy and successful reverse transcription was checked using WTA before continuing with bead reinjection and sorting. Agarose beads containing cellular genomes and transcriptomes were reinjected into droplets to perform single-cell HIV detection PCR. Beads were mixed with PCR reagents to achieve a final concentration of 1× TaqPath Mastermix (Thermo Fisher Scientific, A30866), PEG-6000 (0.5% (w/v)), Tween-20 (0.5% (w/v)), F-127 Pluronic (0.5% (w/v)), BSA (0.1 mg ml⁻¹), HIV gag forward primer (CACTGTGTTTAGCATGGTGTTT, 900 nM), HIV gag reverse primer (TCAGCCCAGAAGTAATACCCATGT, 900 nM) and HIV gag hydrolysis probe (CY5-ATTATCAGAAGGAGCCACCCCACAAGA-3′ Iowa Black RQ, 250 nM)⁶⁸. To generate the final 1× reaction mixture concentration, beads were soaked in 2× PCR master mix on a shaker for 30 min in the dark. Next, the beads were centrifuged and loaded into a 3 ml syringe. The remaining 1× PCR master mix (supernatant) was loaded into a separate 3 ml syringe. Finally, the beads and 1× PCR master mix were reinjected in the microfluidic device to encapsulate the beads into 70 µm droplets⁶⁹. Agarose beads were re-encapsulated in droplets with about 70% loading, which is not accounted for in the detection efficiency calculation. Droplets were collected in 40 µl aliquots in PCR strips and thermocycled as follows: 88 °C for 10 min; then 55 cycles of 88 °C for 30 s and 60 °C for 1 min. After thermocycling, droplets were transferred into a 3 ml syringe for microfluidic sorting.

HIV-DNA⁺ and HIV-DNA⁻ droplets were sorted on the basis of the HIV PCR signal using a concentric sorter as previously described³². For HIV-DNA⁻-sorted samples, we sorted 100 cell equivalents based on the number of genomes per hydrogel bead determined previously, collecting a mixture of HIV-DNA⁻ cell droplets and cell-free droplets. For HIV-DNA⁺-sorted samples, we sorted aliquots of 100 droplets. The sorter was run with the following flow rates: 180 μl h⁻¹ cell droplets, 6,000 μl h⁻¹ bias oil (HFE-7500), 250 μl h⁻¹ spacer oil (HFE-7500) and 3,500 μl h⁻¹ extra spacer oil (HFE-7500). To sort, the 2 M NaCl on-chip electrode was polarized using a high-voltage amplifier at 1,200 V, 4,000 Hz for 15 cycles with 120 μs delay. We sorted into 1.5 ml Eppendorf tubes, removed all but 20 µl of the oil, added 50 µl of distilled nuclease-free water and centrifuged the sample at 20,000g for 5 min, and then stored the samples at −80 °C.

Before performing WTA on sorted HIV-DNA⁺ droplets in each participant, we determined the WTA cycle number that was required to amplify transcriptome cDNA from 100 cells in that participant. Accordingly, we first performed WTA on HIV-DNA⁻-sorted sample aliquots. Sorted HIV-DNA⁻ sample aliquots (frozen at −80 °C) were heated to 60 °C on a heat block for 10 min, mixed carefully by pipet and centrifuged at 20,000g for 5 min. The aqueous layer was then transferred to PCR strips and a WTA PCR reaction was performed using the 1× KAPA HiFi Master mix (Roche, KK2601) and 0.4 μM Smart-seq2 primer (AAGCAGTGGTATCAACGCAGAGT). Sorted material was thermocycled as follows: 95 °C for 3 min; then 18–22 cycles of 98 °C for 15s, 67 °C for 20s and 68 °C for 4 min; then 72 °C for 5 min, with a 4 °C terminal hold. The WTA was performed at three different cycle numbers—18, 20, and 22 cycles. All reactions were subsequently purified using a 1.2:1 ratio of AMPure XP beads (Beckman Coulter), with the final elution performed in 20 µl of nuclease-free water. After WTA, the DNA yield was quantified using the Qubit 4 Fluorometer and DNA size distribution was assayed using a Bioanalyzer 2100 with High Sensitivity DNA chip. On the basis of these results, the HIV-DNA⁺-sorted samples were processed as above using the minimal cycle number required to achieve a concentration of greater than 2 ng µl⁻¹ in 20 µl of elution volume.

Sequencing and read preprocessing

Libraries were prepared from transcriptome material sorted by FIND-seq using the Nextera XT Library Preparation Kit with v2 indexes. Individual sample libraries were combined at equimolar amounts to produce a single library pool. The library was quantified using the KAPA SYBR FAST Universal qPCR Kit. The library concentration and fragment size distribution were confirmed using the Agilent Bioanalyzer 2100 with High Sensitivity DNA chip. The library was diluted and denatured in accordance with the Illumina MiSeq System Denature and Dilute Libraries Guide (document 15039740). Cell line libraries were sequenced on the Illumina MiSeq system in 2 × 75 bp runs, and the selected libraries were subsequently sequenced again on the Illumina HiSeq 4000 system in a 2 × 75 bp run, operated using the Illumina HiSeq Control Software (HCS) v.3.4.0. For samples from people with HIV, libraries were first pooled and run on the Illumina MiSeq system in a 2 × 75 bp run, then rebalanced and run on the Illumina HiSeq 4000 system in a 2 × 75 bp run. Raw sequencing data were converted to fastq format using the bcl2fastq2 script (v.2.20) from Illumina and the reads were demultiplexed using sample-specific indexes. The resulting fastq files were trimmed for quality, ambiguity and presence of read-through adapters using the ‘Trim reads’ tool with the default settings in CLC Genomics Workbench (GWB) v.21.0.3. The quality of the raw and trimmed reads was assessed using QC tools in GWB.

Participant sample data quality filtering

Owing to the abundance of HIV-DNA⁻ cells in samples from ART-treated people with HIV, HIV-DNA⁻ cells were sorted in multiple replicates. Sequencing data were generated from 53 HIV-DNA⁺ and HIV-DNA⁻ cell samples sorted by FIND-seq from 5 people with HIV. A prospective curation approach was used to exclude low-quality samples from downstream transcriptomic analysis. HIV-DNA⁻ sample quality was assessed according to the following parameters: (1) the total number of reads sequenced; (2) the percentage of intergenic and intronic reads; (3) the proportion of ribosomal RNA (rRNA) reads; and (4) the exonic fragment count (Supplementary Table 2). Samples that had a paired-end read count of less than 10⁶ and had >35% mapped intergenic reads were excluded. Furthermore, within each participant, HIV-DNA⁻ samples that differed qualitatively from other replicates by having lower exonic reads or higher rRNA content were removed. If all HIV-DNA⁻ samples were removed for a participant, that participant was excluded from further analysis. After the removal of 31 FIND-seq-sorted samples in this curation process, 22 HIV-DNA⁺ and HIV-DNA⁻ samples belonging to participants 2208, 2510 and 3209 remained (Supplementary Table 2).

Analysis pipeline testing

The transcriptomes of primary cell samples generated by FIND-seq showed high proportions of intronic and intergenic reads (Extended Data Fig. 4). We therefore performed a second, deeper sequencing of libraries from the J-Lat:Raji cell mixing study and tested whether bioinformatics pipelines that address coverage bias and/or genomic DNA contamination might mitigate the effects of these patterns on the gene expression results. In total, we evaluated nine different pipelines using control data from the J-Lat:Raji cell line mixing study. The details of each pipeline are found below; the default options and parameters were used for all tools unless otherwise noted. Reads were mapped against the GRCh38 (ENSEMBL v.100) reference with coding gene annotations only for all pipelines tested.

CLC Genomics Workbench

CLC Genomics Workbench (GWB) v.20 and v.21 (https://digitalinsights.qiagen.com/) were tested using the default settings for mapping and abundance estimation using the RNA-seq analysis tool. For DGE analysis in GWB v.21, the option to filter average expression before FDR correction was selected.

3′ tag counting

Raw reads were preprocessed and mapped using GWB v.21. As in a previous study⁷⁰, reads were mapped to the region within 1,500 bp from the 3′ end of the gene and expression values were calculated in GWB. Analysis of DGE was also performed in GWB.

Salmon with positional bias correction

Salmon v.1.3.0 was implemented as it includes an algorithm for transcript expression quantification that incorporates bias modelling to account for position specific and other biases that are commonly seen in RNA-seq data⁷¹. Read mapping generated from GWB v.20 was used as the input. Post-quantification analysis of DGE was performed using EdgeR (v.3.32.1)⁷² and DESeq2 (v.1.30.1)⁷³.

SeqMonk DNA contamination correction

We considered that relatively high intergenic read proportions in sorted samples might be due to library incorporation of the genomic DNA retained with each cell during FIND-seq. We therefore used the SeqMonk expression quantification (http://www.bioinformatics.babraham.ac.uk/projects/seqmonk/) pipeline v.1.47.2, which estimates and corrects count data for each transcript using the density of intergenic reads. Read mapping previously processed in GWB v.20 was used as the input. Analysis of DGE was performed in DESeq2. Expression qualification and DGE with or without DNA contamination correction (SeqMonk) was evaluated, and each was tested with or without automatic independent filtering (DESeq2).

Selection of the analysis pipeline

For each pipeline, transcriptome accuracy was assessed by comparing J-Lat:Raji FIND-seq mixing study DGE results with the DGE detected between J-Lat cells and the unsorted J-Lat:Raji mixture in standard RNA-seq. DEGs were considered as those with an absolute fold change of ≥1.5 and FDR ≤ 0.05. DEGs identified in standard RNA-seq but not in FIND-seq were considered to be false negatives (FN); those identified only after FIND-seq as false positives (FP); and those identified in both FIND-seq and standard RNA-seq as true positives (TP). Based on this, the sensitivity (or recall) as TP/(TP + FN) and positive predictive value (PPV) as TP/(TP + FP) for each analysis process were calculated (Supplementary Table 7).

GWB v.20 and v.21 yielded the highest combination of sensitivity and PPV. Pipelines that corrected for coverage bias and DNA contamination did not increase the sensitivity, and in several cases showed lower PPV. Although GWB v.20 had a higher PPV than GWB v.21, there were developments in the GWB v.21 transcriptome analysis pipeline that were anticipated to reduce noise in primary cell samples. Thus, the pipeline in GWB v.21 was selected for the analysis of participant samples.

DGE between HIV-DNA⁺ and HIV-DNA⁻ memory CD4 T cells

As described above, transcriptome data from FIND-seq-sorted material contained higher proportions of intronic and intergenic sequences than the standard RNA-seq data. These non-exonic sequences were also abundant in material that was subjected to only the hydrogel encapsulation and cDNA synthesis steps of FIND-seq, consistent with the requisite co-retention of cell genomic DNA with transcriptome material and with efficient nuclear lysis and capture of immature transcripts in our hydrogel-based workflow. Accordingly, after curating the participant samples on the basis of quality, differential expression using only exonic reads was performed (Supplementary Table 3). Using GWB v.21, a combined analysis was performed using the Wald test with Benjamini–Hochberg multiple-testing correction by defining DEGs between HIV-DNA⁺ and HIV-DNA⁻ samples using data from the three participants as biological replicates, while controlling for any interparticipant differences in expression. Moreover, a participant-specific analysis was performed by determining DEGs within each participant separately (Supplementary Table 4). The default settings for all other parameters for the differential expression for RNA-seq tool were used except for Filter on average expression for FDR correction, which was enabled for all analyses. Unless otherwise noted, cut-offs for statistical significance of DEGs were absolute fold change of ≥1.5 and FDR ≤ 0.05.

Euclidean distance calculation

Pairwise Euclidean distances between the curated samples were calculated using the dist function in R (v.4.1.0) using a matrix of counts per million mapped reads (CPM) gene expression values as input. For each sample of a given HIV DNA status group (that is, HIV-DNA⁺ or HIV-DNA⁻), average intragroup and intergroup distances to all other curated samples were calculated, with values plotted in GraphPad Prism (v.9.3.1). Statistical significance of distance differences between groups was calculated using Mann–Whitney U-tests.

Transcriptomic pathway expression differences between HIV-DNA⁺ and HIV-DNA⁻ cells

Ingenuity Pathway Analysis (Qiagen, summer release 2021) was used to identify enriched biological pathways (Supplementary Table 5) on the basis of DEG lists. For the combined analysis considering samples from different participants as biological replicates, DEGs with an absolute fold change of ≥1.5 and FDR ≤ 0.05 were used. For the participant-specific analysis, DEGs with a fold change of ≥2 and raw P ≤ 0.01 were used and pathways regulated in the same direction for all three participants were identified.

The directionality of enrichment of pathways for each analysis was determined from the z-score, which is calculated in Ingenuity Pathway Analysis to represent predicted relative pathway activity. The z-score for each pathway was calculated using the list of genes annotated to that pathway and meeting criteria for statistically significant differential expression between HIV-DNA⁺ and HIV-DNA⁻ cells. A simplified z-score was calculated as follows: Z = (N⁺ − N^–)/(√N), where N⁺ and N^– are those genes of which the direction of regulation is concordant or discordant with predictions from the literature. A positive z-score implies activation of a pathway, whereas a negative z-score implies inhibition. Statistical significance of the enrichment of a pathway was determined using a right-tailed Fisher’s exact test as described previously⁷⁴. Networks of pathways identified as inhibited across participants and their corresponding genes were plotted using ClusterProfiler (v.4.1.1)⁷⁵.

WGCNA

Weighted gene co-expression network analysis⁷⁶ was performed in R using the WGCNA package (v.1.70) with a gene expression matrix of CPM values. Genes detected in <2 samples were excluded from analysis. The one-step automatic method was used for network construction and module detection. A soft thresholding power (β) of 6 was selected based on approximate scale-free topology using the function pickSoftThreshold. The co-expression network was built with a minimum module size of 30, reassignThreshold of 0 and mergeCutHeight of 0.25. The default values were used for the other parameters. Co-expressed modules of genes that correlated with HIV-DNA⁺ and HIV-DNA⁻ status were identified. Modules that were correlated with the traits with P ≤ 0.05 were considered to be significant. GO enrichment analysis for the genes belonging to the two WGCNA modules significantly correlated with cell HIV DNA status was performed using Enrichr (29 March 2021 release)^77,78. Enrichment analysis was performed using a Fisher’s exact test with Benjamini–Hochberg multiple-testing correction.

Analysis of HIV reads

To identify sequence reads representing HIV RNA, we created a combined human (GRCh38, ENSEMBL v.100) and HIV (GenBank: KT284371) reference. The HIV sequence for this reference was derived from the clade B representative in the 2016 LANL HIV sequence compendium, with deletions in the LTR regions replaced by the corresponding sequence and annotations from HXB2CG (GenBank: K03455 M38432), and with masking of the gag amplicon detected in FIND-seq. Reads were aligned to the combined reference using the Map reads to reference tool with the default settings in GWB (v.21). Counts were obtained for reads extracted from mapping to the combined reference. Mapped reads were visualized using GWB and Integrated Genome Viewer (v.2.11.9).

The frequencies of sequence variants in HIV reads compared to the reference sequence were examined to assess the presence of multiple virus sequences. To do this, a consensus of aligned sequences was generated and reads mapping to the HIV genome were extracted. These reads were then mapped against the consensus reference sequence. The resulting mapping was improved by local realignment in areas containing insertions and deletions (indels). Variants were then identified using the ‘low frequency’ variant caller in GWB v.21 with a minimum coverage of 2, minimum count of 1, inclusion of broken reads and without relative read direction filter applied. The default options for the other parameters were used. The list of variants obtained was manually inspected and filtered to remove (1) those with a frequency above 50% (thus representing the predominant sequence rather than a minor variant) and (2) those with read count = 1 or that represented presumptive technical insertions in homopolymeric regions.

Moreover, the Sequences from HIV Easily Reconstructed (SHIVER)⁷⁹ pipeline (v.1.5.8) was tested to create a hybrid reference from de novo assembled contigs of HIV reads for individual samples and closely matched reference sequences. In brief, reads were mapped to the GRCh38 (ENSEMBL v.100) reference using the Map reads to reference tool in GWB v.21 with stringent settings, with the length fraction and similarity fraction parameters set to 0.8. Unmapped reads were then collected and paired reads among them were processed using the de novo assembly tool in GWB (v.21) with the default settings. We also tested the iterative virus assembler (IVA; v.1.0.11) to perform de novo assembly from the unmapped reads using the default settings, but did not recover HIV contigs using this tool. Contig sequences obtained from GWB (v.21) were exported in fasta format and were processed using the SHIVER pipeline with the default settings. A clade B HIV genome obtained from the 2016 LANL sequence compendium was used as a reference.

Enrichment analysis of WGCNA modules in defined CD4 T cell subsets

Viably cryopreserved PBMCs from ART-treated people with HIV were thawed and stained for FACS with LIVE/DEAD Aqua stain (Molecular Probes) and the following antibodies (with the indicated dilutions): CXCR5-Alexa Fluor 488 (1:7; BD), CCR5-Cy7PE (1:10; BD), CD27-Cy5PE (1:10; Beckman Coulter), CD45RO-PE-Texas Red (1:12; Beckman Coulter), CD14-PE (1:80; BD), CD11c-PE (1:40; BD), CD3-H7APC (1:5; BD), CCR7-Alexa Fluor 700 (1:8; BD), CD20-APC (1:5; BD), CD56-APC (1:10; BD), T cell receptor gamma delta (TCR-γδ)-APC (1:5; BD), PD1-Brilliant Violet 711 (1:10; BioLegend), CD8-Qdot 655 (1:200; Invitrogen), CD4-Qdot 605 (1:200; Invitrogen), CD57-Qdot 585 (1:50; Invitrogen) and CCR6-Brilliant Violet 421 (1:10; BD). Stained samples were sorted into CD4 T cell subsets using the FACSAria (BD) system by first gating for single cells that were CD3⁺, Aqua^low and negative for CD11c, CD14, CD20, CD56 and TCR-γδ. The remaining events that were CD4⁺ and CD8⁻ were then collected as naive (CD27⁺CD45RO⁻) or memory CD4 T cell subsets (see memory subset definitions in Extended Data Table 2). Sorted cell subsets were processed for total RNA extraction and whole-transcriptome sequencing as described previously⁶³. The resulting data were processed using the standard pipeline in GWB v.21 using the human reference (GRCh38, ENSEMBL v. 100) with only the coding gene annotations. The resulting CPM values were exported and provided as an input to GSEA (v.4.2.3)^80,81. Enrichment of module 5 and 28 signatures (separated into genes upregulated and downregulated between HIV-DNA⁺ and HIV-DNA⁻ cells) was identified in transcriptome data from each memory CD4 T cell subset (with data from the naive CD4 T cell subset serving as a reference). GSEA was run using the default settings for all of the parameters.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.