Introduction

Each year, approximately two million people are infected with human immunodeficiency virus-1 (HIV-1), and about 1.1 million individuals die from HIV-related illnesses. HIV-1 primarily targets CD4+ T lymphocytes, leading to a gradual depletion of this critical T-cell population and the eventual onset of acquired immunodeficiency syndrome (AIDS) [1, 2]. Although antiretroviral therapy (ART) helps inhibit HIV-1 replication [3], the persistent presence of the virus causes chronic inflammation, contributing to HIV-related morbidity and mortality [4, 5]. Additionally, latent proviruses can be reactivated and causing viral rebound if ART is discontinued. As a result, achieving complete clearance of HIV-1 remains a major goal and public health challenge [6].

During HIV-1 infection, cellular immune responses are disrupted, leading to immune cell dysfunction and alterations in the immune cell landscape. Acute HIV infection (AHI) occurs from 3 weeks to 6 months after transmission, following an initial period when the viral load is undetectable [7]. Once proviruses integrate into the host genome, they impact the recruitment of regulatory factors, thereby affecting gene expression [8–10]. These modifications during AHI play a role in the host cell’s initial immune response and contribute to HIV-1 latency [6, 11, 12]. Studying transcriptomic changes during AHI is crucial to identify primary immune responses and regulatory factors that influence HIV-1’s effect on host immune cell function.

Single-cell RNA sequencing (scRNA-seq) and single-cell Assay for Transposase Accessible Chromatin with high-throughput sequencing (scATAC-seq) are powerful techniques that provide insights into cellular diversity, allowing for high-resolution analysis of gene expression and regulatory landscapes in individual cells [13–17]. Recent advancements in sequencing technologies have driven numerous studies to uncover biological factors contributing to HIV-1 pathogenesis [18–25]. However, challenges remain in studying genetic and epigenetic regulation in HIV-1-infected cells: low HIV-1+ CD4+ T cell frequency in PBMCs limits reliable analysis [19, 26], ex vivo activation in previous studies does not capture in vivo heterogeneity [18, 23, 27]. Additionally, although recent studies have attempted direct ex vivo omics analyses [28–30], integrating these datasets remains difficult due to the varying cell types used in each dataset.

To address these limitations, we conducted single-cell multi-omics analyses on PBMCs from people living with HIV (PLWH) in the early stages of infection. Using a targeted sequencing strategy, we identified and analyzed HIV-1-infected cells at the single-cell level. Based on these findings, we established a gene regulatory network specific to HIV-1-infected cells. This study highlights multiple immune pathways and regulatory factors mediating the effects of acute HIV-1 infection on immune cell function, offering potential therapeutic targets for complete HIV-1 eradication.

Results

Result 1: The single cell transcriptional landscape of HIV-1 infected cells in early stage of infection

To uncover immune mechanisms and transcriptional regulation in HIV-1-infected cells, we performed single-cell multiomics sequencing on PBMCs from nine individuals with early-stage HIV infection (< 6 months) (Supplementary Table 1). Using the 10x Genomics Single-Cell platform, we obtained five scRNA-seq and four snRNA-seq data from single cell multiome datasets. We integrated scRNA-seq and snRNA-seq data, identifying 16 distinct cell populations based on marker genes (Fig. 1A, Supplementary Fig. 1A). Focusing on CD4+ T cells as main targets of HIV-1, we defined four subtypes: CD4 naïve T cells, CD4+ effector memory T cells (CD4 EM cells), T helper (Th)2 cells, and Th17 cells (Fig. 1B, Supplementary Fig. 1B). LDHB and SELL were highly expressed in CD4 naïve T cells, whereas CD4 EM cells were characterized by high expression levels of HLA-DRA and GZMK. Among the activated CD4+ T cells, Th2 cells were identified based on the expression of STAT6 and GATA3 and Th17 cells were annotated based on the expression levels of RORA and STAT3.

Fig. 1.

To identify HIV-1 RNA+ cells, we aligned scRNA-seq reads to the HIV-1 genome. Additionally, we conducted targeted long-read sequencing using the remaining scRNA-seq libraries to improve the sequencing depth of HIV-1 transcripts. We identified 985 HIV-1 RNA+ cells and 3,450 uninfected cells. CD4 naïve and EM cells showed a higher proportion of HIV-1 RNA+ cells than Th2 and Th17 (Fig. 1C, Supplementary Fig. 1C).

Gene expression profiling revealed 334 upregulated and 1,002 downregulated genes in HIV-1 RNA+ CD4+ T cells (Fig. 1D), including genes related to T cell activation (PRDX2, PRDX1, HLA-B, HLA-C, HLA-A, HLA-F, and HLA-E), defense responses to viruses (IFITM3, IFITM1, IFITM2, and MX1), and interferon-beta responses (IFITM3, BST2, IFITM1, IFITM2, XAF1, and SHFL) (Fig. 1E). Notably, IFITM1 was the most upregulated gene in HIV-1 RNA+ CD4+ T cells. As an interferon-stimulated gene, IFITM1 is known to inhibit HIV-1 replication by interfering with viral entry [31]. This finding indicates that the primary cellular pathways active during acute infection are also engaged in the antiviral response to suppress HIV-1.

Given that CD4+ T cell differentiation influences HIV-1 survival [32], we used gene set enrichment analysis (GSEA) to assess subtype-specific responses (Fig. 1F). The GSEA results revealed that interferon-gamma and interferon-alpha responses were upregulated in all CD4+ T-cell subtypes upon HIV-1 infection. Additionally, hallmarks related to T-cell activation were downregulated in both Th2 and Th17 cells. Additionally, while interferon-alpha responses were enriched, downstream pathways (e.g., PI3-AKT-mTOR) were downregulated in Th17 cells (Fig. 1F, G). Downregulated genes in Th17 cells were linked to PI3K-AKT signaling, while no such enrichment was observed among upregulated genes (Supplementary Fig. 2A, 2B). These results suggest that activated T cells experience immune dysfunction in early HIV-1 infection, contrasting with the effects observed in naïve CD4+ T cells.

Result 2: Characterization of epigenetic changes in HIV-1 infected CD4+ T cells in early infection

Single-cell multiome technology enables a detailed understanding of cell states by simultaneously profiling transcriptome and epigenetic signatures within the same cell, facilitating the construction of an accurate gene regulatory network by identifying gene expression drivers. A combined UMAP was generated based on chromatin accessibility and gene expression, with cell types annotated using marker genes (Fig. 2A, Supplementary Fig. 3A). Additionally, coverage plots revealed cell type-specific chromatin accessibility for nearby marker genes (Fig. 2B). Chromatin accessibility in HIV-1 RNA+ naïve and memory CD4+ T cells identified specific patterns, with closed DARs in naïve cells linked to TNF-alpha signaling via NF-κB and open DARs related to HIV-1 Nef protein and signal transduction. In memory CD4+ T cells, closed DARs were related to interferon signaling and negative regulation of viral transcription, whereas open DARs were related to Tat-mediated HIV elongation. When examining genes adjacent to open chromatin regions (Fig. 2D), we observed high chromatin accessibility for PSMB7 and CDK9 in naïve and memory CD4+ T cells, respectively, both of which interact with the HIV-1 Tat protein [33, 34]. Furthermore, the chromatin accessibility of NR4A2 and TRIM62, which negatively regulate HIV-1 [35–37], was reduced in HIV-1 RNA+ cells (Fig. 2E).

Fig. 2.

Result 3: Identification of key transcription factors (TFs) governing HIV-1 infection features

Following HIV-1 infection, host transcription factors (TFs) regulate immune responses and viral persistence [38]. To identify key TFs and construct a gene regulatory network associated with HIV-1 infection, we integrated gene expression and chromatin accessibility data from HIV-1 RNA+ CD4+ T cells.

We found 15 TFs with increased and 24 with decreased regulon activity in HIV-1 RNA+ cells compared to uninfected cells across all CD4+ T-cell subtypes (Fig. 3A). We also confirmed the significantly different regulon activity of these TFs by cell type. CD4 naïve and EM cells shared similar TF activity profiles, while Th17 and Th2 cells exhibited more comparable patterns (Fig. 3B). KLF2, typically high in naïve cells and suppressed in activated cells [39–41], showed elevated regulon activity in naïve and EM cells, with lower activity in Th17 and Th2 cells (Fig. 3B). However, in HIV-1-infected cells, KLF2 activity was upregulated across all subtypes (Fig. 3A). Conversely, FOXO1 showed higher regulon activity in Th17 and Th2 cells but was downregulated in HIV-1 RNA+ Th17 and Th2 cells. In HIV RNA+ cells, the TFs KLF2 and JUND showed strong regulatory activity and elevated gene expression (Fig. 3C). Additionally, we identified the upregulated target genes of six TFs (KLF2, JUND, STAT1, JUNB, XBP1, MYC) in HIV-1 RNA+ cells (Fig. 3D), which were significantly associated with the defense response to viruses, negative regulation of viral genome replication, and innate immune responses (Supplementary Fig. 3B). As the GSEA result showed different enriched signaling pathways in HIV-1 RNA+ cells by various cell types (Fig. 1G), several key transcription factor-encoding genes (KLF2, FOXO1, TBL1XR1) which regulate PI3K-ATK signaling and p53 signaling, were also differentially expressed in each cell types. (Fig. 3E)

Fig. 3.

The chromatin accessibility data further supported the increased activity of these TFs in HIV-1 RNA+ cells. We calculated TF activity based on chromatin accessibility and identified enriched TF-binding motifs in HIV-1 RNA+ CD4+ T cells using DARs. KLF2-, JUND-, and FOSL2::JUNB-binding motifs were highly enriched in the DARs opened in HIV-1 RNA+ cells (Fig. 3F), along with enhanced TF activity (Fig. 3G). Furthermore, TF footprinting analysis revealed higher chromatin accessibility adjacent to the KLF2-binding motif sites (Fig. 3H). The promoter region of KLF2 was more accessible, with higher mRNA expression in HIV-1 RNA+ cells than in uninfected cells (Fig. 3I). Furthermore, the putative cis-regulatory element of KLF2 showed more chromatin-accessible peaks in HIV-1 RNA+ CD4+ T cells than in uninfected cells (Fig. 3I). Overall, this multidimensional approach enabled integrating gene expression and chromatin accessibility in the early stage of infection in PLWH, highlighting potential key TFs that regulate HIV-1 transcription in host cells.

Result 4: KLF2 is an essential regulator in the function of HIV-1-infected CD4+ T cells

Using a multidimensional approach, we identified KLF2 as a key TF within the gene regulatory network of HIV-1 RNA+ CD4+ T cells. KLF2 is reported to be involved in the quiescence, survival, and differentiation of T cells [42]. However, its function in viral infection remains poorly understood. To elucidate the role of KLF2, we curated 88 KLF2 target genes from published chromatin immunoprecipitation (ChIP)-chip, ChIP-seq, and other TF-binding site profiling data, as well as SCENIC regulon datasets. Out of 88 KLF2 target genes, 44 were upregulated in HIV-1 RNA+ cells (Fig. 4A, Supplementary Table 2). These target genes were mainly associated with negative regulation of apoptotic process and negative regulation of programmed cell death (Fig. 4B). Among these significantly upregulated KLF2 target genes, those reported to be involved in the regulation of HIV-1 replication (PIM1 and IL32) [43, 44], cell proliferation (TXNIP) [45], and the JAK/STAT pathway (LTB) also exhibited higher levels of chromatin accessibility in HIV-1 RNA+ CD4+ T cells than in uninfected cells (Fig. 4C).

Fig. 4.

Considering that CD8+ T cells play a critical role in the immune response against HIV-1 infection during chronic infection [46, 47], we further attempted to identify specific ligand– receptor pairs between CD4+ T cells and cytotoxic T lymphocytes to characterize the immune response induced by HIV-1-infected cells [46]. Toward this end, we computed significant ligand–receptor pair interactions between CD8+ and CD4+ T cells using CellPhoneDB (Fig. 4D) [48]. The results indicated a strong interaction between LTB and its receptor (LTBR) in naïve CD4+ T cells, whereas other CD4+ T cell subtypes displayed a weaker interaction with cytotoxic CD8+ T cells. LTB expression was also significantly elevated in HIV-1 RNA+ cells (Supplementary Table 2). Additionally, the NRP2 receptor exhibited heightened levels of interaction between CD4+ and CD8+ T cells, particularly in the effector CD4+ T cell population. NRP2, a known target gene of KLF2, is expressed on T cells and modulates immune responses under pathological conditions [49].

Subsequently, we explored the interactions between the host proteins encoded by KLF2 target genes and HIV-1 proteins. Using the STRING database, we identified 16 KLF2 target genes that encode proteins that directly interact with HIV-1 (Fig. 4E). One of these target genes, SOCS3 which encodes a protein exhibiting a strong interaction with the Tat protein, pivotal for viral genome replication [50]. Another KLF2 target gene, PIM1, which is known as a transactivator of Tat [44], displayed substantial interaction with JUNB and exhibited high gene expression and chromatin accessibility in HIV-1 RNA+ cells (Fig. 4C).

Result 5: Th17 cells exhibit increased enrichment and susceptibility to HIV-1 infection

In vivo analysis of the states of cells with the HIV-1 provirus integrated into host genome has been a long-standing challenge in HIV research. To address this issue, we defined cells with two or more scATACseq reads that aligned to the HIV-1 genome as HIV-1 DNA+ cells. The sequencing reads predominantly aligned to the 5′- and 3′-long terminal repeat (LTR) regions of the HIV-1 genome (Fig. 5A). We observed clustered HIV-1 DNA+ cells in a distinct region, indicating that these cells share a similar epigenetic profile (Fig. 5B, Supplementary Fig. 4A). DARs were identified in the cluster enriched with HIV-1 DNA+ CD4+ T cells compared to the remaining CD4+ T cells (Fig. 5C). The genes closest to the DARs that showed increased accessibility levels in the HIV-1 DNA+ CD4+ T cell cluster compared to the corresponding levels in the uninfected clusters were associated with the inflammatory and defense responses, which was in line with the activation of the innate immune response observed in HIV-1 RNA+ cells (Fig. 5D).

Fig. 5.

Notably, we identified distinctly enriched TF-binding motifs in HIV-1 DNA+ CD4+ T cells and HIV-1 RNA+ CD4+ T cells. The binding motifs for key regulators of Th17 cell differentiation (RORC, STAT3, and RORA) were enriched in HIV-1 DNA+ CD4+ T cells (Fig. 5E, F) [51]. These TFs were also reported to regulate viral gene expression by binding to the HIV-1 LTRs [52]. Overall, HIV-1 DNA+ cells mostly represented the Th17 cell type (Supplementary Fig. 4B, C), in contrast to HIV-1 RNA+ CD4+ T cells, which belonged to various subtypes. Additionally, we confirmed significant enrichment of the KLF2 TF-binding motif in HIV-1 DNA+ CD4+ T cells, as found for HIV-1 RNA+ CD4+ T cells. Furthermore, we detected increased gene expression levels and chromatin accessibility of CCR6, IL4I1, NFKB1, and CCR5 (Fig. 5G).

Discussion

In response to early HIV-1 infection, various host immune responses are triggered, contingent upon the subtype of CD4+ T cells. Our study utilized scRNA-seq and scATAC-seq datasets to profile PBMCs from early HIV-1-infected patients, revealing specific immune responses and mechanisms in distinct HIV-1-infected CD4+ T cell subsets that contribute to HIV pathogenesis. We further identified key regulators of HIV-1 persistence in CD4+ T cells, offering new targets for complete viral elimination.

Interferons play a crucial role in the initial antiviral response during early viral infection by upregulating interferon-stimulated genes and activating signaling pathways like JAK/STAT, PI3K, and p53-dependent pathways. These pathways facilitate T-cell differentiation and selective proliferation of antiviral T cells, effectively inhibiting viral replication and clearing infected cells [53]. However, our findings indicate that during acute HIV-1 infection, Th2 and Th17 cells exhibit heightened interferon-alpha responses, despite diminished downstream signaling pathways. This suggests a potential immunological malfunction induced by HIV-1, hindering host antiviral responses by disrupting interferon target signaling, akin to strategies observed in other viruses like Influenza, HCV, and HSV [54].

In CD4+ memory T cells, we observed enriched chromatin accessibility linked to Tat-mediated HIV-1 elongation, concurrent with downregulated interferon-gamma signaling. This underscores HIV-1’s ability to persist in memory CD4+ T cells and evade host antiviral responses early in infection. Through integrated analysis, we identified KLF2 as a pivotal transcription factor within the regulatory network of HIV-1 RNA+ CD4+ T cells. KLF2, known for promoting cellular quiescence and regulating T-cell migration, exhibited sustained expression across all HIV-1-infected CD4+ T-cell subtypes, including activated cells [39]. Upregulated KLF2 target genes associated with HIV-1 replication and TXNIP, a cell proliferation inhibitor, further suggest HIV-1’s ability to persist without triggering transcriptional repressors or proliferation.

Using scATAC-seq data, we differentiated HIV-1 DNA+ CD4+ T cells, revealing enriched TF motifs (RORC, STAT3, RORA) indicative of Th17 cell predominance. These cells, unlike HIV-1 RNA+ CD4+ T cells distributed across various subtypes, exhibited heightened susceptibility and integration of HIV DNA, with increased expression of CCR6 and CCR5—critical HIV-1 co-receptors [55, 56]. Moreover, enriched KLF2 motifs in HIV-1 DNA+ cells positively regulate CCR5 expression, along with increased IL4I1 expression linked to Th17 cell immunosuppression and NFKB1 involved in HIV-1 genome transcription initiation [57, 58]. This suggests Th17 cells’ potential role as reservoirs due to their susceptibility, quiescent nature, and stem cell-like characteristics.

Limitations include the inability of our HIV-1 targeted sequencing strategy to distinguish HIV-1 RNA from DNA within single cells, necessitating deeper sequencing to capture transcriptionally silenced proviruses. Our study’s design focused on early infection stages (< 6 months), limiting insights into longitudinal changes in cell populations and gene expression profiles during HIV-1 progression. Given the heterogeneity of HIV-1 infection among PLWH, our findings primarily apply to early infection stages and require validation across different infection phases.

In summary, our study unveils disrupted immune mechanisms in distinct HIV-1 RNA+ CD4+ T cell subsets and underscores KLF2’s role in HIV-1-infected cells. We propose that this dysregulation enables HIV-1 transcription maintenance while attenuating pathogenic immune responses. Furthermore, Th17 cells may serve as reservoirs due to their heightened susceptibility and KLF2-mediated CCR5 expression. Our findings provide insights into the complex transcriptional and epigenetic changes during HIV-1 infection, informing potential immunomodulatory therapies to restore normal immune responses and enhance viral control in PLWH, thereby limiting disease progression.

Methods

Preparation of PBMCs samples from early infection of PLWH

PBMCs were collected from a total of 9 individuals in the early stages of HIV infection (<6 months). Epidemiological and virological characteristics of these participants at baseline, including sex, age, duration from HIV-1 infection, CD4 T cell count, HIV-1 viral load and HIV-1 p24 antigen/antibody screening (Supplementary Table 1).

Ethics statement

Written informed consent was given by all participants, and the study adhered to the principles of the Declaration of Helsinki and Korea Disease Control and Prevention Agency (IRB No. 2019-07-06-P-A) approved this study.

Preparation of scRNA-seq libraries

The libraries were prepared using the Chromium Single Cell 3′ Library & Gel Bead Kit v2 (PN-120237) (10x Genomics) following the manufacturer’s instructions. The single-cell suspension was washed twice with 1x PBS containing 0.04% BSA. Cell number and viability were determined using the Bio-Rad TC20 cell counter. The cells were loaded onto the 10x Genomics Chromium Controller to generate gel beads in emulsions (GEMs). Library preparation was performed according to the 10× Genomics Chromium Single Cell 3′ reagent kit (V2 chemistry) instructions. The quality and concentration of the libraries were assessed using the Agilent Bioanalyzer 2100. Finally, the libraries were sequenced on an Illumina NovaSeq.

Preparation of scATAC-seq libraries

The libraries were prepared using the Chromium Next GEM Single Cell ATAC Library & Gel Bead Kit (PN-1000175) following the manufacturer’s instructions. After isolating the nuclei from PBMCs, the single-nuclei suspension was mixed with Transposition Mix and incubated for 60 minutes at 37°C. Subsequently, the transposed nuclei were combined with the master mix and loaded onto the 10x Genomics Chromium Controller to generate gel beads in emulsions (GEMs). Library construction was performed according to the instructions provided with the Chromium Next GEM Single Cell ATAC Reagent Kits (v1.1). The subsequent steps were performed in the same way as described in ‘Preparation of scRNA-seq libraries’.

Preparation of single cell multiome libraries

The libraries were prepared using the Chromium Next GEM Single Cell Multiome ATAC + Gene Expression Reagent Kits (PN-1000283) (10x Genomics) following the manufacturer’s instructions. Briefly, nuclei were mixed with Transposition Mix and incubated for 60 minutes at 37°C. Then these were loaded onto Chromium Next GEM Chip J (10x Genomics) and GEMs were generated according to manufacturer’s instructions. After post-GEM cleaned up, barcoded transposed DNA and barcoded full length cDNA fragments were pre-amplified with PCR. ATAC libraries were constructed using the Library Construction Kit (PN-1000190) and full-length pre-amplified cDNA fragments were additionally amplified via PCR. RNA libraries were constructed using the Library Construction Kit B (PN-1000279). Finally, the libraries were sequenced on an Illumina NovaSeq 6000.

HIV-1 targeted sequencing

After scRNA-seq, remaining scRNA-seq full-length cDNA libraries were subjected to HIV-1 targeted sequencing using the SQK-LSK-109 Ligation Sequencing Kit (Oxford Nanopore Technologies) and three specific primers for HIV-1 (Supplementary Table 3). The TrueSeq Read1 primer was used to preserve the cell barcode information, and the cDNA fragments containing the HIV-1 sequence were amplified by the HIV-1 specific primers. The cDNA fragments with HIV-1 were primarily selected through PCR using biotinylated primers. Subsequently, the HIV-1 specific primers selectively amplified the cDNA fragments through hemi-nested-PCR. These selected cDNA fragments were sequenced on a MinION FLO-MIN106 R9.4.1 flow cell (Oxford Nanopore Technologies). Base calling was performed by using Guppy (V6.1.1), and cells with HIV-1 reads were identified using the cell barcode information.

Sequence alignment to human and HIV-1 genome

The reads from scRNA-seq and scATAC-seq were aligned to GRCh38 human genome (v.3.0.0) using CellRanger count (10x Genomics, v.6.1.2) and CellRanger ARC (v.2.0.0) with default options, respectively. Feature-barcode matrices were used for downstream analysis. To identify the sequencing reads containing the HIV-1 sequence, the reads from scRNA-seq were aligned to the HIV-1 genome using the STARsolo (v. 2.7.1a). We adjusted the mapping percentage parameter from 50 % to 25% to map the short length of HIV-1 sequence included in the sequencing reads.

scRNA-seq data processing

For the scRNA-seq datasets, Seurat R package (v.3.2.1) was used for processing the raw count matrices (UMI counts per gene per cell).[59] Cells with fewer than 200 features and greater than 30% mitochondrial gene expression were excluded, and cells expressing fewer than three genes were eliminated for downstream analysis. Subsequently, the dataset was normalized, highly variable genes were selected, and scaling was performed in Seurat. Batch effects were removed using the ‘RunFastMNN’ function, and principle component analysis (PCA) was performed after integrating the datasets. Cells were clustered based on the computed nearest neighbor map using the ‘FindNeighbors’ function. The clustered cells were visualized on UMAP embeddings. Differentially expressed genes in each cluster were identified using the ‘FindAllMarkers’ function, and the clusters were annotated using the known cell markers. Clusters were merged if the number of differentially expressed genes was less than 10.

scATAC-seq data processing

For the scATAC-seq datasets, Signac R package (v.1.5.0) was used to process the data.[60] To ensure common features across the datasets, combined peaks were generated. The transcription start site (TSS) enrichment score and nucleosome signal score were calculated for each cell using the ‘TSSEnrichment’ and ‘NucleosomeSignal’ functions in Signac. Cells were retained with a TSS enrichment score > 1 and a nucleosome signal < 2, and further filtered the cells containing < 1000 or > 100,000 ATAC fragments. MACS2 was utilized to identified peaks using each fragment file of datasets.[61] The datasets were merged, and batch effects were removed using Harmony.[62] The merged dataset was subjected to dimensional reduction with LSI using the “RunTFIDF”, “FindTopFeatures” and “RunSVD” functions. The LSI components were used for UMAP projection. The cell type of each cluster was annotated using differentially accessible peaks. To identify significant peaks for each cluster, the ‘FindAllMarks()’ function was employed with the following parameters: min.pct = 0.1, test.use = ‘LR’, only.pos = TRUE. The ‘ClosestFeature()’ function was then used to find the closest genes to each peak within each cluster. Cell clusters were annotated using known cell markers. Subsequently, the scRNA-seq dataset and scATAC dataset from the 10x multiome dataset were integrated using identical cell barcodes to ensure the scATAC-seq clusters’ cell type annotation, leveraging the previously annotated scRNA-seq cell types.

Gene ontology enrichment analysis and gene set enrichment analysis

To identify enriched biological pathways in HIV-1 RNA+ CD4 T cells, DEGs for each cell type were identified using ‘FindMarkers’ function in Seurat (v.3.2.1). We compared the normalized gene expression data of HIV-1 RNA+ CD4 T cells with uninfected CD4 T cells. DEGs were filtered based on |logFC| > 0.25 and adjusted P-value < 0.01. Gene ontology enrichment analysis was performed using the web-based tools of DAVID (version DAVID 6.8; https://david.ncifcrf.gov/). Biological processes (BPs) were ranked based on the sets of DEGs in HIV-1 RNA+ cell, considering BPs with a P-value < 0.01 significant. Gene set enrichment analysis (GSEA) was performed to identify the statistically significant gene sets in the HIV-1 RNA+ cell group within different CD4 T cell subtypes. The raw count matrix was used for the analysis. Hallmark gene sets from Molecular Signatures Database (MSigDB) were used. Only gene sets containing 10 to 500 genes were selected.

Gene regulatory network analysis

To construct a gene regulatory network and predict TF activities from a gene expression dataset, SCENIC (R package, https://github.com/aertslab/SCENIC) was utilized as a computational method.[63] Co-expression modules between TFs and target genes were identified using GENIE3. The cisTarget Human motif database was utilized for scoring the regulons. Subsequently, the AUCell method was employed to score single cells and identify enrichment of regulons in the infection group.

Cell-cell interaction analysis

CellphoneDB was utilized to analyze cell-cell communication between CD4 T cell subtypes and CD8 T cell subtypes.[48] Interaction pairs with a positive mean expression level in the HIV-1 RNA+ group were selected to explore associations among cell types. The interaction strength was visualized using the ‘plot_cpdb’ function in the ktplots R package, represented as a dot plot.

Transcription factor binding motif enrichment and footprinting analysis

The cells were categorized into three groups: HIV-1 RNA+ cell group, HIV-1 DNA+ cell group, and uninfected cell group. Differential accessibility regions (DARs) among the groups were identified using the ‘FindMarkers’ function. DARs with a p-value less than 0.01 between the cell groups were filtered. DARs with a positive log2 fold change value were selected for further analysis to identify enriched motifs. Enriched motifs were identified using the motifmatchr (R package, https://bioconductor.org/packages/motifmatchr), which involved testing for overrepresentation of each DNA sequence motif in the DARs using a hypergeometric test. The computation of GC contents and matching of background sets of peaks for GC content, sequence length, and counts were performed using the ‘FindMotifs’ function in Signac. The per-cell motif activity score was calculated by chromVAR for each motif across the infection conditions.[64] For the identification of a strong enrichment of Tn5 integration events adjacent to the enriched motifs, TF footprinting analysis was conducted using the ‘Footprint’ function in Signac.

Statistical analysis

The subsection above provides a description of statistical analyses for single cell transcriptomic and epigenetic studies. All descriptive statistical analyses were conducted in R 3.6.0 and higher. Benjamini-Hochberg correction was used for significance.

Data availability statement

The raw data was deposited in Korean Nucleotide Archive (KoNA, https://kobic.re.kr/kona)* with the accession ID, KAP230707.

Acknowledgements

We thank all study participants by donating their blood and made this study possible. We are grateful to Jaehyun Seong for preparing this study samples of acute HIV-infected patients. The study samples of acute HIV-infected patients are collected by the Chronic Infectious Disease Cohort Study (Grant Number 4800-4859-304). Also this work was supported by the Korea National Institute of Health (Grant Number: 2019-NI-067 and 2020-ER5103).

Additional information

Contributors

Conceptualization: JP, BSC

Data collection: SYC, SSI

Data curation: DL

Formal analysis: DL

Investigation: CHY, SIS, HA

Funding acquisition: JP, BSC

Supervision: JP, BSC

Writing – original draft: DL

Writing – review & editing: JP, BSC, DL

Funding

This work was supported by the Korea National Institute of Health (Grant Number: 2019-NI-067 and 2020-ER5103).

Additional files

Supplementary material.