Integrated Single-Cell Analysis of Multicellular Immune Dynamics during Hyper-Acute HIV-1 Infection

Cellular immunity is critical for controlling intracellular pathogens, but the dynamics and cooperativity of the evolving host response to infection are not well defined. Here, we apply single-cell RNA-sequencing to longitudinally profile pre- and immediately post-HIV infection peripheral immune responses of multiple cell types in four untreated individuals. Onset of viremia induces a strong transcriptional interferon response integrated across most cell types, with subsequent pro-inflammatory T cell differentiation, monocyte MHC-II upregulation, and cytolytic killing. With longitudinal sampling, we nominate key intra- and extracellular drivers that induce these programs, and assign their multi-cellular targets, temporal ordering, and duration in acute infection. Two individuals studied developed spontaneous viral control, associated with initial elevated frequencies of proliferating cytotoxic cells, inclusive of a previously unappreciated proliferating natural killer (NK) cell subset. Our study presents a unified framework for characterizing immune evolution during a persistent human viral infection at single-cell resolution, and highlights programs that may drive response coordination and influence clinical trajectory.

acute infection, we performed scRNA-Seq on peripheral blood mononuclear cells (PBMCs) from 96 four individuals enrolled in FRESH who became infected with HIV, assessing multiple timepoints 97 from pre-infection through one year following initial detection of viremia (Fig. 1A, table S1). In our 98 study, hyper-acute infection refers to timepoints at and prior to peak-viral load, whereas acute 99 infection refers to timepoints after peak viral load but before 6 months. Samples were processed 100 in duplicate using Seq-Well 18 -a portable, low-input massively-parallel scRNA-Seq platform

108
To assign cellular identity, we performed variable gene selection, dimensionality reduction, 109 clustering, and embedding en masse across data collected from all individuals and timepoints 110 (see Methods). Samples were combined for cell type/phenotype identification to find common 111 transcriptional features of ubiquitous cell subsets, and to improve statistical power on classifying 112 small/rare cell types. Importantly, combined analyses yielded few individual-specific features in 113 the resulting clustering and embedding, suggesting that disease biology, rather than technical 114 batch, is the main driver of variation and subsequent clustering (Fig. 1D, Fig. S1A,B). We 132 average CD4 + -R 2 = 0.491, p = 0.0416; average CD8 + -R 2 = 0.665, p = 0.00158 (Fig. 1E).

133
Subsequently, we calculated frequencies for the other cell types in our scRNA-Seq data as a

152
To identify tightly co-regulated modules (M) of genes for each type for each individual, we 153 applied weighted gene correlation network analysis (WGCNA) 26,27 on all cells classified as a 154 particular cell type across all timepoints ( Fig. 2A; see Methods for details). Strongly correlated 155 gene modules (permutation test for within-module similarity, FDR corrected q < 0.05) were then 156 tested for significant variation over time by scoring cells at each timepoint against the genes within 157 a module, followed by tests for shifts in score distribution between pairs of timepoints (Wilcoxon 158 rank sum test, FDR corrected q < 0.05). This generated 0-8 temporal modules per cell type (for a 6 list of all significant modules see Table S3 for gene membership and Table S4 for

175
Beginning with one individual (P1), we identified a set of 6 significant gene modules spanning 176 multiple cell types that all shared their highest relative module score at the peak viremia timepoint 177 (Fig. 2B). Upon inspection of the genes within each, we uncovered a core set of genes shared

233
We subsequently grouped these MMs by temporal shape (    Monocytes expressed genes associated with antigen presentation and IL-4 signaling (mainly 257 HLA-DR subunits), which may reflect generalized interferon responses, or the potential to 258 promote active T helper and CTL responses. NK cells, CTLs, and proliferating T cells all 259 upregulated genes associated with killing of target cells by perforin and granzyme release, 260 highlighting the joint role of innate and adaptive cells in combating viremia (see Table S5

266
To identify common and cell-type specific inducers of these measured transient responses 267 extending past peak viremia, we generated a list of predicted upstream drivers of each module 268 (see Table S6). Selecting highly significant hits in at least two modules, we drew a network of

282
We also contextualized observed responses to these upstream drivers temporally by re-283 scoring against enriched genes for each driver. This analysis revealed variable kinetics in the 284 onset, intensity, and length of immune responses across different cell types (Fig. 3G, Fig. S7).

285
We note the following gene-programming upregulation trends in most individuals: CD4+ T cells

305
After discovering temporally variant modules in our dataset, we observed a few sets that 306 demonstrated similar temporal response patterns in a given cell type, but were not combined into 307 a single module by our framework. We thus sought to understand how these modules might be 308 linked by looking across single cells for module co-expression. Here, single-cell expression data 309 are essential to distinguish response circuitry among cells.

310
The clearest example of multiple modules being co-expressed with the same temporal 311 pattern in the same cell type from our analysis was the NK activated M3 module (CCL3, CCL4, 312 CD38) and the cytotoxic M4 module (PRF1, GZMB, HLA-A) in P3 (Fig. 3D), both part of MM1.

313
Enrichment analysis demonstrated little overlap between the significant pathways associated with 314 these modules, implying orthogonal biological function. We therefore investigated whether the 315 gene programs for these modules were highly co-expressed in the same single cells and thus co-316 varied among single cells across time (Fig. S8A). While we did not observe differential 317 simultaneous upregulation of these modules between time points, we found variation in the 318 correlation of cell-scores for the pair as a function of time across single cells, with the strongest 319 correlation one to two weeks after HIV detection (Fig. S8B). Variation in the correlation of M3 and 320 M4 may reflect a synergizing of these gene programs 77 within NK cells to combat HIV as viremia 321 declines post peak.

11
In examining MM3 (Fig. S5) -containing the majority of the IFN response modules -we 323 observed that P3 also exhibited a set of temporally similar modules in monocytes (M1 & M3); 324 however, these modules did not variably correlate in expression score as a function of time.

325
Instead, these gene programs were highly co-expressed but only at HIV-detection (Fig. S8C-D).

326
Gene set analysis readily demonstrated that monocyte M1 consisted of IFN response genes, 327 while M3 was enriched for genes associated with inflammation (Fig. S8E). IFN has been shown 328 to stunt the production of pro-inflammatory cytokines in monocytes similar to the phenotype 329 observed in these cells in viremic persons 78,79 , but the co-expression of anti-viral and pro-330 inflammatory signals in the same single cells has not yet been described to our knowledge. As

349
In all four individuals, we saw dramatic structuring of the monocytes in PC space by time.

350
Specifically, monocytes sampled at HIV detection (0 weeks) or 1-week post-detection were 351 strongly separated along either PC1 or PC2, indicating a pervasive hyper-acute response to 352 infection. Interestingly, non-classical monocytes (see Fig. S1D and Table S2), which may be more 353 susceptible to infection 80 , displayed disparate temporal dynamics across individuals, even though 354 they drove significant variation in PCA space (Fig. S9D). Comparing DE genes at these peak 12 response timepoints (vs. pre-infection) highlighted not only the specificity of the co-356 inflammatory/anti-viral monocytes to P3, but also other person specific differences in monocyte 357 phenotype (Fig. 4C). Gene set analysis on upregulated genes in each individual confirmed that    Table S7).

392
Recently, we demonstrated that, in most individuals in the FRESH study, a majority of 393 proliferating CTLs in hyper-acute infection are HIV-specific by tetramer staining 86 . Therefore, we 394 turned to the proliferating T cells in our study to look for differences in response based on long-395 term viral control. En masse, the proliferating T cells expressed similar levels of cytotoxic genes 396 as non-proliferating CTLs (Fig. S10C). DE analysis highlighted genes associated with cell-cycle 397 (e.g. STMN1, HIST1H1B, MKI67) and memory (e.g. IL7R, KLRB1) (see Fig. S10D

412
We next utilized unsupervised analyses to explore differences in proliferating T cell 413 responses over time among individuals (Fig. 5B, Fig. S10F). Proliferating T cells captured at 1-

420
To determine whether these TRDC + CD16 + cells were gdT or NK cells, we scored them, as well as 14 non-proliferating CTLs and NK cells, against gene signatures described in that study (Fig. S10G).

422
Based on score similarity to NK cells, and the relative down-regulation of CD3 compared to the