Single Cell Analysis of Blood Mononuclear Cells Stimulated Through Either LPS or Anti-CD3 and Anti-CD28

Immune cell activation assays have been widely used for immune monitoring and for understanding disease mechanisms. However, these assays are typically limited in scope. A holistic study of circulating immune cell responses to different activators is lacking. Here we developed a cost-effective high-throughput multiplexed single-cell RNA-seq combined with epitope tagging (CITE-seq) to determine how classic activators of T cells (anti-CD3 coupled with anti-CD28) or monocytes (LPS) alter the cell composition and transcriptional profiles of peripheral blood mononuclear cells (PBMCs) from healthy human donors. Anti-CD3/CD28 treatment activated all classes of lymphocytes either directly (T cells) or indirectly (B and NK cells) but reduced monocyte numbers. Activated T and NK cells expressed senescence and effector molecules, whereas activated B cells transcriptionally resembled autoimmune disease- or age-associated B cells (e.g., CD11c, T-bet). In contrast, LPS specifically targeted monocytes and induced two main states: early activation characterized by the expression of chemoattractants and a later pro-inflammatory state characterized by expression of effector molecules. These data provide a foundation for future immune activation studies with single cell technologies (https://czi-pbmc-cite-seq.jax.org/).


INTRODUCTION
The immune system plays a central role not only in fighting infections, but also in the pathogenesis of many diseases. While access to disease-relevant tissues may be limited, blood profiling provides a minimally invasive method to assess immune function and health (1,2). Blood serves as a conduit for transport of immune cells throughout the body and circulating immune cells frequently display signatures of disease or immune dysfunction (i.e., immunosenescence) (3)(4)(5). Transcriptional profiling of blood provides an unbiased method for immune monitoring and has proved invaluable in uncovering novel disease mechanisms. In vitro activation of immune cells GRAPHICAL ABSTRACT | Overview of the study design and analysis methods.
is an effective way to assess their functional capacity. For example, T cells can be activated with non-specific mitogens such as PHA (phytohemagglutinin) and ConA (Concanavalin A) (6), or more specifically through the T cell receptor and co-stimulatory receptors, either by crosslinking with antibodies against CD3 and CD28 (7), or with antigen-laden antigen presenting cells (8). Similarly, innate cells, including monocytes, can be activated through pattern recognition receptors (PRRs) with pathogen-associated molecular patterns (PAMPs) such as lipopolysaccharide (LPS) (9). Studying such cellular responses is an effective strategy to uncover functional defects and/or disease-associated phenotypes in human immune cells, and has been widely applied to autoimmune diseases (10)(11)(12), immunodeficiency, aging (13)(14)(15)(16), and allergy (17).
Advances in single cell technologies can precisely uncover donor-, cell type-, and cell state-specific variation in immune cell functions from clinical samples containing mixed cell populations (18,19). Hence, these technologies are increasingly used to study ex vivo and in vitro cellular responses of diverse immune cell types [CD8 + T cells (20,21), dendritic cells (DCs) (22,23), and B cells (24)] as well as activated cells in immune diseases (25). However, single cell gene expression profiling of activated cells has unique data generation and analysis challenges including identifying and removing multiplet cells that increase in number as cellular interactions increase upon activation. Furthermore, avoiding batch effects is essential to distinguish activation-specific gene expression patterns from batch-induced changes. Another major limitation is the high cost of the methodology which prevents large scale studies. Recent advances can address some of these challenges. For example, multiplexing via hashtag oligonucleotides (26) and individuals' genotypes (27) permits simultaneous sequencing of libraries from multiple individuals and conditions to reduce both the cost and technical batch effects. Finally, cellular indexing of transcriptomes and epitopes by sequencing (CITEseq) technology (26) profiles epitopes (i.e., cell surface proteins) and gene expression levels from the same single cells and is particularly suitable for studying distinct immune cell types and states (e.g., activated and resting) and for filtering out multiplets among activated cells.
We combined cell hashing and genotype-based multiplexing to generate CITE-seq data for the deep analysis of circulating human immune cell responses to classic activators of adaptive (anti-CD3/CD28) and innate (LPS) immune cells. Responses to these activators have been studied over decades (6) via bulk profiling of sorted cells in a reductionist manner, which missed indirect effects of these activators on other cell types as well as the heterogeneity in responses at the single cell level. PBMCs from 10 healthy young donors were profiled at baseline and following activation via (1) anti-CD3 + anti-CD28 for selective T cell activation, or (2) LPS for monocyte activation using CITEseq to quantify changes in both transcript and protein levels using 39 antibodies (e.g., CD3, CD4, CD14, CD56, CD69, CD25). Data from this study uncovered cell-specific responses to each condition at both the gene and protein level and uncovered direct and indirect effects of these activators on immune cell transcriptomes. This study will guide future single cell studies of activated immune cells in health and disease (https://czi-pbmccite-seq.jax.org/).

Simultaneous Protein and Gene Expression Profiling of Resting and Activated PBMCs
PBMCs were isolated from the blood of 10 healthy young individuals (5 men, 5 women; 21-32 years old; Caucasian), and were cultured for 24 h with LPS (10 ng/ml), anti-CD3/CD28 (2 µg/ml each) (section Methods), or in medium alone (referred to as "Baseline"). Flow cytometry analysis confirmed that anti-CD3/CD28 treatment induced the expression of activation markers (CD25, CD69) on T cells (Figures 1A,B), whereas LPS activated monocytes as shown by CD80 induction (Figures 1C,D). To study the effects of these activation conditions at the single cell level, we pooled baseline and activated cells from the 10 individuals (30 samples in total) and used CITE-seq to simultaneously profile changes in levels of transcripts and 39 cell surface proteins (Supplementary Table 1). To reduce experimental batch effects and costs, equal numbers of cells from all donors and conditions were pooled and sequenced together. Cells from different conditions were labeled with Cell Hashing Antibodies (BioLegend) (Figure 1E), where three hashtag oligonucleotides (HTOs) were used to demultiplex data from each experimental condition (baseline, LPS, anti-CD3/CD28) using a Gaussian mixture model on HTO counts (Supplementary Figure 1A). Finally, donor genotypes were used to demultiplex data from each individual with Demuxlet (27) (Figure 1F) Overall, 195,480 droplets were captured from all samples (Supplementary Table 2). Cell multiplets from different sources were detected and eliminated with the help of multiplexing strategies (HTOs and genotypes) and the levels of cell-specific marker proteins (Supplementary Figure 1B, section Methods). HTOs uncovered multiplets from different treatment conditions, which constituted ∼13% (26,084/195,480) of total droplets in addition to 18,871 (∼9%) empty droplets (i.e., no HTO detected) (Supplementary Figure 1C). Using genotypes, we identified and filtered out multiplets from different individuals as well as cells that could not be assigned to an individual with high confidence (i.e., ambiguous) (Supplementary Figure 1D), resulting in 55,353 single cells with a mean of 1,932 unique molecular identifiers (UMIs) per cell. Empty droplets from cell hashing and ambiguous cells from Demuxlet overlapped extensively (83%, Supplementary Figure 1E). Finally, expression levels of cell type-specific proteins uncovered multiplets from different cell types (section Methods). After multiplet removal and demultiplexing, an average of 1,638 (range 914-2,034) cells per donor were detected (16,382 total cells) ( Figure 1G), with 855 genes detected per cell on average. Donors showed similar cell frequency distributions across conditions ( Figure 1H). The proportions of CD25 + /CD69 + T cells from CITE-seq and flow cytometry data were highly correlated (Pearson's R = 0.96) (Supplementary Figure 1F) across individuals upon anti-CD3/CD28 stimulation, indicating that CITE-seq is an effective technology to capture and study immune cell activation.

Anti-CD3/CD28 Treatment Activates All Classes of Lymphocytes
Baseline and anti-CD3/CD28-activated PBMCs were jointly analyzed using the top 500 most variably expressed genes and proteins (Figure 2A; section Methods). Cell clusters were annotated using the expression of marker proteins: CD19 and CD20 for B cells, CD16 and CD56 for NK cells, CD14 and CD11c for CD14 + monocytes, CD3 and CD4 for CD4 + T cells, CD3 and CD8a for CD8 + T cells, and CD45RA and CD45RO for naïve and memory T cells, respectively ( Figure 2B and Supplementary Figure 2A). Accordingly, baseline PBMCs from 10 individuals consisted of 60-82% T cells, 12-25% CD14 + monocytes, 4-15% NK cells, and 3-6% B cells (Figure 2C), which was highly correlated with flow cytometry data from this cohort (Supplementary Figure 2B; Pearson's R = 0.97). CD8 + memory T cells were the most variable cell type between individuals with a coefficient of variation (CV) of 0.65, followed by NK cells (CV = 0.51).
Anti-CD3/CD28 treatment activated all lymphocytes (Figure 2A, left). PBMCs grouped primarily by cell type and secondarily by cell state (i.e., resting vs. activated) based on their annotations with cell surface markers ( Figure 2B and Supplementary Figure 2A). Labeled antibodies allowed the detection of naive (CD45RA high ) and memory (CD45RO high ) T cells. Activated T cells showed early (CD69) and late (CD25) activation markers and reduced CD3 expression ( Figure 2B). Indirect activation of NK and B cells was revealed by the induction of CD69 mRNA and protein ( Figure 2B). Notably, this treatment caused a significant decline in monocyte numbers and their percentages in PBMCs ( Figure 2C). This is likely due to the induction of apoptosis in monocytes via activated T cells (28,29) as apoptosis-associated genes were upregulated (e.g., CTSL, TRAF1, TNF, and GZMB) among T cell-monocyte multiplets detected in the anti-CD3/CD28 condition.
Differential expression analyses between activated and resting T cells revealed that 950 genes were upregulated in naïve and memory CD4 + and CD8 + T cells (Supplementary Figure 2D), including IFNG, Type I interferon-inducible genes (ISG15, ISG20), and interferon regulatory factors (IRF1, IRF4) (Supplementary Table 4). Stimulation induced more genes in naïve (CD45RA + CD45RO − ) T cells than memory (CD45RO + CD45RA − ) T cells; 824 transcripts were induced specifically in naïve CD4 + and CD8 + T cells, while 494 were induced in CD8 + naive T cells (Supplementary Figure 2D). Naïve cell-specific response genes were enriched for KEGG pathways associated with cell cycle (e.g., MCM2/3/5/6, CDKN2D, SMC1A, CHEK1), carbon metabolism (e.g., ESD, CS, SDHD, GPI), and basal transcription factors (e.g., ETF1, EIF4E, EIF4H, EIF4G1), emphasizing the proliferative and expansive capacity of naive T cells upon activation (Supplementary Table 4). Pseudo-temporal ordering of baseline and activated T cells using the Slingshot algorithm (34) revealed mostly linear trajectories for activation despite the existence of distinct memory (CD45RO + CD45RA − ) subsets at baseline, thereby indicating that these different subsets respond similarly to anti-CD3/CD28 activation. Distinct subsets (e.g., CD161 + CCR6 + ) were no longer detectable upon activation-potentially due to internalization of cell surface proteins-preventing further interrogation of their specific responses. Instead, our analyses revealed that, for both CD4 + T and CD8 + T compartments, activated naïve (CD45RA + CD45RO − CD69 + ) cells were ordered toward the end of the activation pseudotime trajectories, due to the greater transcriptomic changes in naive T cells upon anti-CD3/CD28 stimulation compared to memory cells (Figures 3D,E). Taken together these data show that anti-CD3/CD28 activation has a significant effect on T cell transcriptomes, where naïve T cells go through more significant changes compared to their memory counterparts.

Anti-CD3/CD28 Activation of PBMCs Result in Indirect Activation of NK Cells
Anti-CD3/CD28 treatment induced CD69 protein expression on NK cells, and transcription of 74 genes (Supplementary Table 5), 60 of which were shared with T cells, including cytotoxic molecules (e.g., GZMB) as well as interferon response genes (Supplementary Figure 2D, Supplementary Table 4). Clustering analyses of NK and T cells, at baseline and after anti-CD3/CD28 activation, uncovered the transcriptional similarity between NK and CD8 + T cells ( Figure 4A). Interestingly, activated CD8 + T cells clustered distinctly from NK cells likely due to the upregulation of distinct molecules in T cells (Supplementary Figure 2D). To further characterize these cells, we derived MAIT, cytotoxicity, and senescence scores using relevant gene sets (section Methods, Supplementary Table 3). As expected, CD8 + CD161 + CCR6 + cells had the highest MAIT score. Activated NK cells (CD56 + CD69 + ) showed the highest cytotoxicity score among all studied cell types, followed by activated CD8 + memory cells (CD8 + CD45RO + CD69 + ) T cells ( Figure 4B). In alignment with this, upon activation, important cytotoxic molecules were up-regulated in NK cells, including perforin (PRF1), granulysin (GNLY), and granzymes (GZMB) (Supplementary Table 5). In contrast, activated T cells had high senescence scores, particularly for activated naive CD8 + T cells (CD45RA + CD69 + ), which was in alignment with the increased expression of CDKN1A (encoding p21 protein) and CDKN2A (encoding p16 protein) along with other senescence-associated molecules ( Figure 3C). In contrast, baseline and activated NK cells displayed low senescence scores.
The significant effect of anti-CD3/CD28 treatment on the NK cell transcriptome (Figure 4C; left) results from indirect activation, since neither TCR nor CD28 are expressed by NK cells. Unsupervised clustering of NK cells from baseline and anti-CD3/CD28 stimulation identified seven groups (Figure 4C; middle), four of which were enriched in activated cells (∼95% activated, Clusters 0, 1, 4, and 5) (Figure 4C; right). Cells in these four clusters of activated NK cells had higher expression of interferon-stimulated genes (ISGs) including LY6E, IFI6, and IFI27 ( Figure 4D) than NK cells at baseline. One of these clusters (cluster 4) highly expressed marker proteins for activation: CD69 and CD83 (35) (Figure 4D). NK cells in cluster 4 expressed XCL1 that facilitates their communication with DCs, as well as CCL3 and CCL4 (36) that enable them to recruit various immune cells including naive CD8 + T cells. In contrast, clusters 2, 3, and 6 cells were mostly from the baseline condition and enriched for KLF2, which encodes a transcription factor important for NK cell survival (37). Trajectory analyses of NK cells (baseline and activated) revealed that the cells at the end of the trajectory have the highest expression of interferon-gamma (IFNG), granzyme B (GZMB), CD69, and various chemokines (CCL3, CCL4, XCL1, XCL2) (Figure 4E), resembling the profiles of cluster 4 cells. Expression of other interferon-inducible genes (IFITM1, ISG15) and granulysin (GNLY) peaked prior to the end of the trajectory, highlighting differences in timing and degree of NK subtype activation. Thus, scRNAseq of anti-CD3/CD28 activated PBMCs revealed that NK cells are indirectly activated with this condition that resulted in upregulation of cytotoxic molecules.  Table 5). Increased transcription included ribosome and proteasome activity, oxidative phosphorylation machinery, and mitochondrial stress pathways. It also involved genes associated to humoral immune responses and Type 2 interferon signaling, such as CD40, CXCL9, CXCL10, and IRF1 (Supplementary Table 5).

Anti-CD3/CD28 Activation of PBMCs Results in B Cell Activation
Resting and activated B cells clustered into six distinct groups ( Figure 5B; middle), where clusters 0 (comprising 25% of B cells) and 3 (comprising 15% of B cells) were specifically enriched in activated B cells (>90%) (Figure 5B; right). Cells from the activated B cell clusters highly expressed ISGs (IF44, IF44L, STAT1, and MX1) and the activation marker CD69 (Figure 5C). Cluster 0 cells expressed genes associated with age-associated B cells (ABC) (38)(39)(40) and DN2 cells (41) (e.g., TLR7, PRDM1, TGAX/CD11c and TBX21/T-bet). Distinct from ABCs and DN2 cells, they retained follicular (CXCR5) and conventional memory (CD27) markers. Cluster 3 cells expressed additional markers associated with DN2 cells, such as FCRL5 and TRAF5, but lacked the expression of TLR7, ITGAX, and TBX21. Differences detected between these activated B cells and ABC/DN2 cells may reflect acute vs. chronic activation of B cells as seen with aging and auto-immune diseases. Slingshot analyses revealed mostly linear activation trajectories for B cells, where the cells at the end of the trajectory had the highest expression of activation markers (CD25, CD69, and CD278) as well as genes encoding ribosomal proteins ( Figure 5D). Thus, scRNAseq of anti-CD3/CD28 activated PBMCs is providing new clues to the generation of DN2 cells, that are expanded in SLE (25,41).

LPS Treatment of PBMCs Only Affects CD14 + Monocytes and Induces Inflammation-Associated Genes
Baseline and LPS-activated PBMCs were jointly analyzed using the top 500 most variably expressed genes and proteins ( Figure 6A) and were annotated using the expression of marker proteins (section Methods) (Supplementary Figure 3A). LPS stimulation did not lead to significant changes in PBMC cell compositions ( Figure 6B). Among PBMCs, only monocytes were activated by the LPS treatment (Figures 6A,C) and showed upregulated expression of pro-inflammatory cytokine genes (e.g., IL8, IL1B) (Supplementary Figure 3B). To globally assess the enrichment of inflammation-associated transcripts in each cell type, we used a reference list of 249 inflammationrelated genes from NanoString Technologies Inc. (section Methods). Activated monocytes had significantly higher inflammation scores compared to baseline monocytes and all other cell types ( Figure 6D). Overall, 219 genes were differentially expressed in monocytes (125 up-and 94 downregulated; FDR 5%) upon LPS treatment (Figure 6E and Supplementary Table 5). Functional enrichment analyses of these differentially expressed genes (DEGs) using gene ontology (GO) terms and KEGG pathways (42) revealed that LPS treatment activated pattern recognition receptor signaling pathways (Toll-like and NOD-like) and the proinflammatory NF-kappa B signaling pathway (FDR 5%; Figure 3C), including the up-regulation of important chemokines (CCL3, CCL4) and cytokines (IL1B, IL6) ( Figure 6F). In contrast, LPS stimulation decreased the expression of genes associated with antigen processing and presentation via MHC class II and phagosome activity (Figure 6F, Supplementary Table 6).

LPS-Activated Monocytes Exhibit Two Distinct Transcriptional States
The Slingshot (34) algorithm revealed that the trajectory of monocytes followed a linear transition from resting to activated states (Figure 6G; top, Supplementary Figure 3D). LPS-activated monocytes showed increased variation in their transcriptional responses and were positioned at the end of the trajectory (Figure 6G; top). Genes that are important for defining monocyte trajectories included inflammatory cytokines (IL6, IL8, IL1B), metallothionein genes involved in zinc transport [MT1F (43)] and oxidative stress [MT2A (44)], as well as alarmins (S100A8/S100A9), which induce pro-inflammatory cytokine production (45) (Supplementary Figure 3E). Alarmins and metallothionein gene expression peaked in early stages (pseudotime values 5 up to 15), whereas pro-inflammatory cytokines (e.g., IL1B, IL6, IL8) peaked in late stages (pseudotime ≥ 15) ( Figure 6G; bottom, Supplementary Figure 3E). These LPS-induced cellular states and changes were observed consistently across donors (Supplementary Figure 3F). To further dissect this heterogeneity, we clustered baseline and activated monocytes using k-means clustering (k = 4). Baseline monocytes clustered separately from LPS-treated ones ( Figure 7A, Supplementary Figure 4A), which were split into 3 clusters referred to as LPS1, LPS2, and LPS3 states. These clusters coincided with the pseudotime values (Supplementary Figure 4B), where LPS1 and LPS2 states mapped to the early activation state and LPS3 mapped to the later inflammatory state (Supplementary Figure 4B). These four monocyte clusters were detected in all individuals ( Figure 7B); inflammation scores of clusters (section Methods) increased in agreement with their position in the trajectory, LPS3 having the highest inflammation score (Figure 7C). LPS2 and LPS3 states had the most distinct transcriptional response profiles, with 50 Frontiers in Immunology | www.frontiersin.org

Inflammation scores of all cells under two conditions. LPS-activated monocytes have the highest expression levels of inflammation-associated transcripts. (E)
Differential expression analysis of LPS-stimulated and baseline monocytes revealed 125 induced and 94 reduced genes, respectively (FDR 5%), as well as 3 induced and 3 reduced proteins, respectively. (F) Enriched KEGG pathways (FDR 5%) for LPS response genes. Node size for each pathway represents the number of genes in the pathway that were also induced in LPS-stimulated monocytes. Node color for genes indicate induction (red) or reduction (blue) upon LPS treatment. LPS-activated monocytes upregulate genes associated with pattern recognition receptor signaling pathways (Toll-like receptor, NOD-like receptor) and inflammation (NF-kB signaling), but also downregulate genes associated with antigen presentation and phagosome activity. (G) Trajectory (pseudotime) inference of monocytes under baseline and LPS conditions identified genes with heterogeneous activity in stimulated monocytes. Cells are color-coded with respect to activation status (top left) and inferred pseudotimes (top right). Expression levels of selected genes were overlaid on cells sorted in the activation trajectory (bottom). and 47 genes uniquely induced in these states, respectively, compared to baseline (Supplementary Figures 4C,D,  Supplementary Table 7).
Network analyses of genes induced in LPS2 and LPS3 states using Ingenuity Pathway Analysis (IPA; section Methods) revealed that LPS2 genes included alarmins (S100A8, S100A9) and genes associated with cell movement and trafficking ( Figure 7D; genes in bold text, pink outline), such as various chemoattractants CXCL7 (PPBP), CCL7, and CCL8. In contrast, LPS3 genes included pro-inflammatory interleukins (e.g., IL1A, IL1B, IL6), macrophage inflammatory protein-1-alpha (CCL3),−1-beta (CCL4), and−3 (CCL20), and tumor necrosisassociated factors (TNFAIP6, TNIP1) ( Figure 7E). Interestingly, in LPS3, we also observed the induction of anti-inflammatory genes including IL10, HLA-E-an MHC class I molecule with inhibitory effects on cytokine production in monocytes (46), IDO1-an immunosuppressive enzyme -and one of its regulators PTGS2 (47), suggesting a feedback loop to control increased inflammation. Pathway enrichment analysis using IPA disease and function terms verified that pro-inflammatory signals emerge at LPS2 state and peak at the later LPS3 state (Supplementary Figure 4E). At baseline, monocytes showed higher relative expression of ribosomal protein genes (Supplementary Figure 4F) and genes involved in mRNA processing, translation, and protein synthesis. Thus, scRNAseq analysis further delineates how LPS activates monocytes into early and late inflammatory states.

DISCUSSION
In vitro activation of circulating immune cells have been widely used for immune monitoring and identification of disease mechanisms since 1980s (8,48). However, such assays are generally limited in scope and often only provide information regarding a single cell type. A holistic study of circulating immune cell responses to activators will uncover how the system as a whole respond to an immune challenge and how cell-cell interactions contribute to these responses. Hence these studies will bring us closer to a better understanding of immune cell responses in health and disease. CITE-seq enables precise studies of in vivo and in vitro immune cell responses. However, these assays are costly to generate and harbor technical challenges that might confound data interpretation (i.e., multiplets, batch effects). Herein we describe a cost-effective and in-depth study to uncover responses of peripheral blood mononuclear cells (PBMCs) from 10 healthy donors (5 men, 5 women) to two widely used activators-anti-CD3/CD28 and LPS-to activate adaptive and innate cells, respectively. PBMCs were profiled using CITE-seq technology with 39 antibodies before and after activation; libraries were multiplexed using both cell hashing and genotypes, which permitted the super-loading of cells to reduce costs and batch effects (49). This design enabled studying activated cells in a cost-effective manner while also detecting and labeling multiplets from (i) different conditions, (ii) different individuals; and (iii) different cell types. We observed an increase in biological multiplets upon activation, potentially due to increased cell-cell interactions, which should be taken into consideration in future study designs to decouple biological (i.e., cellular interactions) from technical (i.e., due to super-loading) multiplets (50,51). One caveat of multiplexing many libraries is the limited number of cells detected per donor (∼1,600 cells), which was enough to study responses of major cell populations (e.g., CD14 + monocytes), but not that of less abundant ones (e.g., CD16 + monocytes, dendritic cells). Multiplexing fewer samples can effectively address this caveat in future studies.
Anti-CD3/CD28 treatment activated all classes of lymphocytes either directly (T cells) or indirectly (B, NK cells). Upon activation, co-stimulatory molecules (CD28, CD278) and activation markers (CD69, CD25) were expressed on the surface of both naïve and memory T cells. Transcriptomic changes in naïve T cells were the greatest likely due the fact that these cells start from a "lower transcriptional state" compared to their memory counterparts as they need to differentiate and initiate clonal expansion (52). Activated T cells highly expressed senescence (e.g., CDKN1A-encoding p21, CDKN2A-encoding p16, SESN2, and PRKAA1), effector and cytotoxic molecules (e.g., GZMA/B/H, IFNG). Using CITE-seq antibodies, under baseline conditions, we detected distinct subsets of memory cells (e.g., Th subsets, MAIT). However, these subsets were not detected upon activation, potentially due to the internalization of cell surface molecules and increased transcriptional similarity between memory subsets.
Indirect activation of NK cells via anti-CD3/CD28 activated 74 genes, most of which were shared with T cell responses including CD69, interferon stimulated/response, and cytotoxic genes. Scoring of single cells using different gene sets showed that activated NK cells had the highest scores for cytotoxicity whereas activated T cells had highest scores for senescence. Clustering and trajectory analyses revealed heterogeneity among activated NK and B cells. Late stage activated NK cells resembled previously described "helper" NK cells and expressed CD83 protein (35), whereas late stage activated B cells expressed molecules that have been expressed in B cells associated with aging and autoimmune diseases (TLR7, PRDM1, ITGAX/CD11c, TBX21/T-bet) (39)(40)(41). However, compared to DN2s and ABCs, activated B cells retained follicular and conventional memory markers (CXCR5 and CD27), likely as a reflection of acute vs. chronic activation. Future single cell activation studies from different populations should help further characterize these B cells. LPS stimulation specifically activated monocytes and pseudotemporal analyses uncovered two states of activated monocytes: (i) a first state that includes pro-inflammatory alarmins (S100A8/A9) and molecules to attract and interact with other immune cells (PPBP, CCL7, CCL8); and (ii) a second, possibly more terminated inflammatory state, with increased expression of bona fide inflammatory cytokines (IL1B, IL6, IL8, and IL24) and potential anti-inflammatory molecules (IL-10, SOCS3, HLA-E, and IDO1) to mitigate the effects of increased inflammation. A similar phenomenon was noted in human blood dendritic cell subsets in which cellular responses to influenza virus infection resulted in three waves of chemokine secretion to attract effector immune cells (53). Induction of interferon response genes and TFs (STAT1 and IRF1) previously reported in IFN-g-primed (24 h) (54) and LPS-induced (0-4 h) monocytes (43), was not observed in our study, potentially due to the 24-h stimulation, which is a later time-point in the activation cascade. Anti-CD3/CD28 treatment reduced the numbers of monocytes in PBMCs likely due to the induced apoptosis via activated T cells (28). LPS treatment had no effect on B cells because B cells from healthy human donors lack TLR4 on their surfaces (55)(56)(57)(58)(59) and respond poorly to TLR4 (60,61).
In summary, this study provides an experimental and computational framework to carefully study activated immune cells using single cell methods by multiplexing to handle batch effects and cut costs and by quantifying protein expression to identify cell types, cell states, and multiplet cells. The data, analyses, and application developed here will be important resources for future studies such as exposure of immune cells to antigens, microbes, adjuvants and other targeting agents. It will guide the design of subsequent experiments to further expand our understanding of cellular responses to additional in vivo and in vitro activation conditions in healthy and diseased individuals. Overall, this study provides a foundation for future immune activation experiments using single cell technologies. The data is freely shared via Human Cell Atlas platform as well as via an interactive R shiny application (https://czi-pbmc-cite-seq.jax. org/).

Blood Collection
Blood was collected at the UConn Health Clinical Research Center (CRC) at 263 Farmington Avenue Farmington, CT 06030-3805 following informed consent of participants for a UConn Health IRB-approved protocol (IRB Number: 18-151J-2). Approximately 50 ml of blood was collected from each donor. Samples were collected BD vacutainer tubes containing Acid Citric Dextrose (ACD) Solution A, and peripheral blood mononuclear cells (PBMC) were isolated immediately after.

Isolation and Stimulation of PBMCs
Peripheral blood mononuclear cells (PBMCs) were isolated with gradient centrifugation using Lymphoprep (StemCell Technologies) and the red blood cells were removed by incubating the cells with RBC lysis buffer (Tonbo Biosciences). After that, the PBMCs were immediately cultured using activation conditions described below for 24 h at 1 million cells/ml density in complete RPMI (RPMI supplemented with 10% heat-inactivated Fetal bovine serum, 100 U/ml Penicillin/Streptomycin, 2 mM L-glutamine, 10 mM HEPES, 0.1 mM Non-essential amino acids, 1 mM Sodium pyruvate) at 37 • C + 5% CO 2 . The PBMCs were cultured in three different ways: (1) Plate-bound anti-CD3 + anti-CD28 (clones OKT3 and CD28.2 at 2 µg/ml each-Tonbo Biosciences)-plates were coated with Goat anti-mouse IgG (clone Poly4053 at 10 µg/ml-BioLegend) for 1 h at 37 • C, washed, then coated with both anti-CD3 and anti-CD28 also for 1 h at 37 • C. After that cells were added to the plates for stimulation; (2) lipopolysaccharide (LPS) (10 ng/ml-InvivoGen); and (3) Complete medium only as control. After 24 h in culture, an aliquot of each sample was stained for flow cytometry analysis and the remaining cells were harvested and frozen for subsequent CITE-seq analysis and genotyping.

Sample Demultiplexing Using Cell Hashing
Three hashtag oligonucleotides (HTO) were used to multiplex samples from different experimental conditions (e.g., HTO A for baseline samples; HTO B for Anti-CD3/CD28 treated samples). Cell HTOs were counted in raw sequencing data using CITEseq-Count version 1.4.2. (49). CITE-seq-Count was used to count HTOs for each cell using recommended parameter settings ("cells 25,000 -cbf 1 -cbl 16 -umif 17 -umil 26 -max-error 3"). HTO counts per cell were normalized with respect to the total read counts for that cell and scaled and log2 transformed. Distributions of read counts for each HTO across all cells were bimodal. To classify whether a cell has detectable counts for an HTO, we fit a Gaussian mixture model (with 2 mixture components) for each HTO using mclust R package (64) (version 5.4.5). Based on these models, we labeled each cell as: (i) empty: if there are no sufficient read counts for any HTO; (ii) singlet: only if counts for a single HTO were detected; (iii) multiplets: if there are detectable counts for 2 ("doublet") or 3 ("triplet").

Sample Demultiplexing With Demuxlet
Genotypes were used to demultiplex the donor of origin for each sample using Demuxlet version 2 (27) with default parameters. For this we used 167,703 genotyped SNPs overlapping exonic regions in the genome. Cells annotated to one individual were retained, while "ambiguous" and "doublet" cells were discarded from further analyses.

RNA Quantification
Single cell RNA-sequencing libraries were prepared using the Chromium Single Cell 3 ′ Reagent Kit (v2 Chemistry) according to the user guide (10× Genomics). Resulting libraries were sequenced using Illumina NovaSeq S4 flowcells to an average depth of ∼150,000 reads per cell. Reads were aligned to human reference sequence GRCh37/hg19 and unique molecular identifiers (UMIs) were quantified using the Cell Ranger "count" software (10x Genomics, version 2.1.0) and with additional parameters "-expect-cells=25000." Cell barcodes with fewer than 400 detected genes, and those not identified to be singlets via cell hashing or Demuxlet were discarded from analyses.

Protein Quantification
Read counts for antibody derived tags (ADTs) for 39 distinct antibodies were measured from raw sequencing reads using CITE-seq-Count with the following parameters: "-cells 25,000cbf 1 -cbl 16 -umif 17 -umil 26 -max-error 3". Raw counts were normalized using a centered log ratio approach using Seurat R package version 3.0.2 (65).

Multiplet Cell Removal Using Surface Protein Expression
Number of multiplet cells increases upon activation, due to the increase in cell-cell interactions. To remove such multiplets, we took advantage of cell surface protein expression levels and identified cells that co-expressed multiple cell-type-specific proteins. First, we determined which cells were "positive" and "negative" for each protein using the HTODemux function from Seurat (65). Since some immune cells co-express certain cell surface markers (e.g., NK and CD8 + T cells both express CD8), we devised a multiplet removal protocol for each cell type. For B cell multiplets, we filtered out cells that express CD19 in addition to one of the following (i) a T cell marker (CD3, CD4, CD8); (ii) an innate cell marker (CD14, CD16, CD56). For CD14 + monocyte multiplets, we removed cells that express CD14 and (i) a T cell marker (CD3, CD8, CD28); (ii) NK marker CD56; (iii) a B cell marker (CD19, CD20). For NK cell multiplets, we removed cells that are positive for CD56 and one or more of the following (i) T cell markers (CD3, CD4, CD28); (ii) monocyte markers (CD14, CD11c); (iii) B cell markers (CD19, CD20). For CD4 + T cell multiplets, we removed cells that express CD4 markers (CD3, CD4) and one or more of the following (i) CD8 (for CD8 + T), CD56 (for NK), CD19, CD20 (for B cells), CD14, CD16 (for monocytes/NK). Lastly, for CD8 + T cell multiplets, we removed cells that co-expressed CD8 markers (CD3, CD8) and one or more of the following (i) CD4 (for CD4 + T cells), (ii) CD56 (for NK), (iii) CD19, CD20 (for B cells), and iv) CD14, CD16 (for monocytes/NK). After all cells are removed we obtained comparable number of single cells for each individual upon 3 conditions. These cell numbers correlated significantly with cell numbers inferred from flow cytometry data (R = 0.97).

Cell Clustering
Cells were clustered using top 500 most variable genes and proteins, which are determined by fitting a line to the relationship of the variance and mean using a local polynomial regression (loess). Protein/gene expression values were standardized using the observed normalized mean and expected variance (given by the fitted line). Feature variance was then calculated on the standardized values after clipping to a maximum of the square root of the number of cells. After scaling and centering the feature values, principal component analysis (PCA) was performed on the top 500 variable features. Following this, uniform manifold approximation and projection (UMAP) was applied on the top 20 principal components to reduce the data to two dimensions. Next, a shared nearest neighbor (SNN) graph was created by calculating the neighborhood (Jaccard index) overlap between each cell and its 20 (default) nearest neighbors. To identify cell clusters, we used a SNN modularity optimization based clustering algorithm (66).

Differential Expression
To identify cell-type-specific response genes and proteins, we conducted pairwise Wilcoxon Rank Sum tests between each PBMC cell-type under stimulation and baseline conditions (e.g., LPS stimulated monocytes vs. baseline monocytes). In each differential expression analysis, genes detected in at least 10% of each cell group were considered. Features with an absolute difference in expression > log2(0.25) and false discovery rate <5% were regarded as significant results.

Single Cell Pseudo-Temporal Trajectory Inference
Benchmarking and inference of single cell pseudotime trajectories were performed with the R package dyno (version 0.1.1) (67). Briefly, for each cell-specific response (e.g., monocytes induced by LPS) trajectories were inferred using Slingshot (the method determined to minimize error and maximize scalability and accuracy). Pseudo-temporal ordering of cells were derived using the corresponding baseline and activated cells as well as the genes and ADTs previously determined to be differentially induced (FDR < 5%) between the two groups of cells.

Ingenuity Pathway Analysis
Response genes were uploaded to Qiagen's Ingenuity Pathway Analysis (IPA) version 01-16 system to build gene networks for each gene set. Networks were constructed using the Connect tool with "Human" as the species parameter and otherwise default parameters to map direct and indirect interactions between molecules as described in the Ingenuity Knowledge Base (IKB). IPA's core analysis function was also used to annotate each response gene set and resultant networks with canonical pathways, relevant diseases, and biological functions.

Single Cell Scoring
Genes associated with inflammation were obtained from NanoString Technologies, Inc. (cite, https://www.nanostring. com/products/gene-expression-panels/gene-expression-panelsoverview/ncounter-inflammation-panels) and used to calculate an inflammation score for PBMCs at baseline and LPSstimulated conditions. Feature counts for each cell were divided by the total counts for that cell, multiplied by a scale factor of 10,000, and then log normalized. Normalized read counts for all inflammation-associated genes were divided by the normalized counts of all transcripts and then min-max scaled between 0 and 1. Gene lists (Supplementary Table 3) were used to score cytotoxicity, MAIT or senescence expressions across T and NK cell clusters. To do so, we calculated the mean expression for each cell, within each cluster.

Pathway Analysis of Stimulation-Response Genes
Differentially expressed genes for each cell-type-response to stimulation were functionally annotated using the R package clusterProfiler (version 3.14.0) (42). The Kyoto Encyclopedia of Genes and Genomes (KEGG), gene ontology (GO), and WikiPathways databases were used to determine association with particular biological processes, diseases, and molecular functions. The top pathways with an FDR-adjusted p < 5% were summarized in the results.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by UConn Health IRB-approved protocol (IRB Number: 18-151J-2). The patients/participants provided their written informed consent to participate in this study.