Single cell analysis of blood mononuclear cells stimulated through CD3 and CD28 shows collateral activation of B and NK cells and demise of monocytes

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/). Graphical abstract


31
The immune system plays a central role not only in fighting infections, but also in the pathogenesis of many diseases. While 32 access to disease-relevant tissues may be limited, blood profiling provides a minimally invasive method to assess immune 33 function and health (1,2). Blood serves as a conduit for transport of immune cells throughout the body and circulating 34 immune cells frequently display signatures of disease or immune dysfunction (i.e., immunosenescence) (3-5). 35 Transcriptional profiling of blood provides an unbiased method for immune monitoring and has proved invaluable in 36 uncovering novel disease mechanisms. In vitro activation of immune cells is an effective way to assess their functional 37 capacity. For example, T cells can be activated with non-specific mitogens such as PHA (phytohemagglutinin) and ConA 38 (Concanavalin A) (6), or more specifically through the T cell receptor and co-stimulatory receptors, either by crosslinking 39 with antibodies against CD3 and CD28 (7), or with antigen-laden antigen presenting cells (8). Similarly, innate cells, Differential expression analyses between activated and resting T cells revealed that 950 genes were upregulated in naïve 157 and memory CD4 + and CD8 + T cells (Fig. S2D), including IFNG, Type I interferon-inducible genes (ISG15, ISG20), and 158 interferon regulatory factors (IRF1, IRF4) (Table S4). Stimulation induced more genes in naïve (CD45RA + CD45RO -) T 159 cells than memory (CD45RO + CD45RA -) T cells; 824 transcripts were induced specifically in naïve CD4 + and CD8 + T cells, 160 while 494 were induced in CD8 + naive T cells (Fig. S2D). Naïve cell-specific response genes were enriched for KEGG 161 pathways associated with cell cycle (e.g., MCM2/3/5/6, CDKN2D, SMC1A, CHEK1), carbon metabolism (e.g., ESD, CS, 162 SDHD, GPI), and basal transcription factors (e.g., ETF1, EIF4E, EIF4H, EIF4G1), emphasizing the proliferative and 163 expansive capacity of naive T cells upon activation (Table S4). Pseudo-temporal ordering of baseline and activated T cells 164 using the Slingshot algorithm (34) revealed mostly linear trajectories for activation despite the existence of distinct memory 165 (CD45RO + CD45RA -) subsets at baseline, thereby indicating that these different subsets respond similarly to anti-166 CD3/CD28 activation. Distinct subsets (e.g., CD161 + CCR6 + ) were no longer detectable upon activation -potentially due to 167 internalization of cell surface proteins -preventing further interrogation of their specific responses. Instead, our analyses 168 revealed that, for both CD4 + T and CD8 + T compartments, activated naïve (CD45RA + CD45RO -CD69 + ) cells were ordered 169 towards the end of the activation pseudotime trajectories, due to the greater transcriptomic changes in naive T cells upon 170 anti-CD3/CD28 stimulation compared to memory cells (Figs. 3D,E). Taken together these data show that anti-CD3/CD28 171 activation has a significant effect on T cell transcriptomes, where naïve T cells go through more significant changes 172 compared to their memory counterparts. 173 174

monocytes and induces inflammation-associated genes 224
Baseline and LPS-activated PBMCs were jointly analyzed using the top 500 most variably expressed genes and proteins 225 ( Fig. 6A) and were annotated using the expression of marker proteins (Methods) (Fig. S3A). LPS stimulation did not lead 226 to significant changes in PBMC cell compositions (Fig. 6B). Among PBMCs, only monocytes were activated by the LPS 227 treatment (Fig. 6A, C) and showed upregulated expression of pro-inflammatory cytokine genes (e.g., IL8, IL1B) (Fig. S3B). 228 To globally assess the enrichment of inflammation-associated transcripts in each cell type, we used a reference list of 249 229 inflammation-related genes from NanoString Technologies Inc. (Methods). Activated monocytes had significantly higher 230 inflammation scores compared to baseline monocytes and all other cell types (Fig. 6D). Overall, 219 genes were 231 differentially expressed in monocytes (125 up-and 94 down-regulated; FDR 5%) upon LPS treatment ( Fig. 6E; Table S5). 232 Functional enrichment analyses of these differentially expressed genes (DEGs) using gene ontology (GO) terms and KEGG 233 pathways(42) revealed that LPS treatment activated pattern recognition receptor signaling pathways (Toll-like and NOD-234 like) and the pro-inflammatory NF-kappa B signaling pathway (FDR 5%; Table S6, Fig. S3C), including the up-regulation 235 of important chemokines (CCL3, CCL4) and cytokines (IL1B, IL6) (Fig. 6F). In contrast, LPS stimulation decreased the 236 expression of genes associated with antigen processing and presentation via MHC class II and phagosome activity (Fig. 6F, 237 Table S6). 238

Discussion 269
In vitro activation of circulating immune cells have been widely used for immune monitoring and identification of disease 270 mechanisms since 1980s (8,48). However, such assays are generally limited in scope and often only provide information 271 regarding a single cell type. A holistic study of circulating immune cell responses to activators will uncover how the system 272 as a whole respond to an immune challenge and how cell-cell interactions contribute to these responses. Hence these studies 273 will bring us closer to a better understanding of immune cell responses in health and disease. CITE-seq enables precise 274 studies of in-vivo and in-vitro immune cell responses. However, these assays are costly to generate and harbor technical 275 challenges that might confound data interpretation (i.e., multiplets, batch effects). Herein we describe a cost-effective and 276 in-depth study to uncover responses of peripheral blood mononuclear cells (PBMCs) from 10 healthy donors (5 men, 5 277 women) to two widely used activators -anti-CD3/CD28 and LPS -to activate adaptive and innate cells, respectively. PBMCs 278 were profiled using CITE-seq technology with 39 antibodies before and after activation; libraries were multiplexed using 279 both cell hashing and genotypes, which reduced costs and batch effects. This design enabled studying activated cells in a 280 cost-effective manner while also detecting and labeling multiplets from i) different conditions, ii) different individuals; and 281 iii) different cell types. We observed an increase in biological multiplets upon activation, potentially due to increased cell-282 cell interactions, which should be taken into consideration in future study designs to decouple biological (i.e., cellular 283 interactions) from technical (i.e., due to superloading) multiplets (49,50). One caveat of multiplexing many libraries is the 284 limited number of cells detected per donor (~1600 cells), which was enough to study responses of major cell populations 285 (e.g., CD14 + monocytes), but not that of less abundant ones (e.g., CD16 + monocytes, dendritic cells). Multiplexing fewer 286 samples can effectively address this caveat in future studies. 287 288 Anti-CD3/CD28 treatment activated all classes of lymphocytes either directly (T cells) or indirectly (B, NK cells). Upon 289 activation, co-stimulatory molecules (CD28, CD278) and activation markers (CD69, CD25) were expressed on the surface 290 of both naïve and memory T cells. Transcriptomic changes in naïve T cells were the greatest likely due the fact that these 291 cells start from a 'lower transcriptional state' compared to their memory counterparts as they need to differentiate and 292 initiate clonal expansion (51). Activated T cells highly expressed senescence (e.g., CDKN1A -encoding p21, CDKN2A -293 encoding p16, SESN2, and PRKAA1), effector and cytotoxic molecules (e.g., GZMA/B/H, IFNG). Using CITE-seq 294 antibodies, under baseline conditions, we detected distinct subsets of memory cells (e.g., Th subsets, MAIT). However, 295 these subsets were not detected upon activation, potentially due to the internalization of cell surface molecules and increased 296 transcriptional similarity between memory subsets. 297 Indirect activation of NK cells via anti-CD3/CD28 activated 74 genes, most of which were shared with T cell responses 299 including CD69, interferon stimulated/response, and cytotoxic genes. Scoring of single cells using different gene sets 300 showed that activated NK cells had the highest scores for cytotoxicity whereas activated T cells had highest scores for 301 senescence. Clustering and trajectory analyses revealed heterogeneity among activated NK and B cells. Late stage activated 302 NK cells resembled previously described 'helper' NK cells and expressed CD83 protein (35) In summary, this study provides an experimental and computational framework to carefully study activated immune cells 320 using single cell methods by multiplexing to handle batch effects and cut costs and by quantifying protein expression to 321 identify cell types, cell states, and multiplet cells. The data, analyses, and application developed here will be important 322 resources for future studies such as exposure of immune cells to antigens, microbes, adjuvants and other targeting agents. 323 It will guide the design of subsequent experiments to further expand our understanding of cellular responses to additional 324 in vivo and in vitro activation conditions in healthy and diseased individuals. Overall this study provides a foundation for 325 future immune activation experiments using single cell technologies. The data is freely shared via Human Cell Atlas 326 platform as well as via an interactive R shiny application (https://czi-pbmc-cite-seq.jax.org/). 327 328 Methods 329

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

Isolation and stimulation of PBMCs 336
Peripheral blood mononuclear cells (PBMCs) were isolated with gradient centrifugation using Lymphoprep (StemCell 337

Technologies) and the red blood cells were removed by inclubating the cells with RBC lysis buffer (Tonbo Biosciences). 338
Atfer that, the PBMCs were immediately cultured using activation conditions described below for 24 hours at 1 million 339 cells/ml density in complete RPMI (RPMI supplemented with 10% heat-inactivated Fetal bovine serum, 100U/ml 340 Penicillin/Streptomycin, 2mM L-glutamine, 10mM HEPES, 0.1mM Non-essential amino acids, 1mM Sodium pyruvate) at 341 37°C + 5% CO2. The PBMCs were cultured in three different ways: 1) Plate-bound anti-CD3 + anti-CD28 (clones OKT3 342 and CD28.2 at 2μg/ml each -Tonbo Biosciences) -plates were coated with Goat anti-mouse IgG (clone Poly4053 at 10μg/ml 343 -BioLegend) for 1 hour at 37°C, washed, then coated with both anti-CD3 and anti-CD28 also for 1 hour at 37°C. After that 344 cells were added to the plates for stimulation; 2) lipopolysaccharide (LPS) (10ng/ml -InvivoGen); and 3) Complete medium 345 only as control. After 24 hours in culture, an aliquot of each sample was stained for flow cytometry analysis and the 346 remaining cells were harvested and frozen for subsequent CITE-seq analysis and genotyping. manufacturer's instructions. We mapped the Illumina array probe sequences to hg19 genome assembly and excluded 369 potential problematic ones as in (56,57). Furthermore, we filtered out SNPs with an alternate allele frequency difference 370 with 1000G EUR samples > 20%, or palindromic SNPs with a minor allele frequency > 20%, genotype missingness > 2.5%, 371 Hardy-Weinberg p-value < 10-4. Finally, 1,147,545 genotyped SNPs were considered for further analyses. an average depth of ~150,000 reads per cell. Reads were aligned to human reference sequence GRCh37/hg19 and unique 393 molecular identifiers (UMIs) were quantified using the Cell Ranger "count" software (10X Genomics, version 2.1.0) and 394 with additional parameters "--expect-cells=25000". Cell barcodes with fewer than 400 detected genes, and those not 395 identified to be singlets via cell hashing or Demuxlet were discarded from analyses. 396 397

Multiplet cell removal using surface protein expression 403
Number of multiplet cells increases upon activation, due to the increase in cell-cell interactions. To remove such multiplets, 404 we took advantage of cell surface protein expression levels and identified cells that co-expressed multiple cell-type-specific 405 proteins. First, we determined which cells were "positive" and "negative" for each protein using the HTODemux function 406 from Seurat (59). Since some immune cells co-express certain cell surface markers (e.g., NK and CD8 + T cells both express 407 CD8), we devised a multiplet removal protocol for each cell type. For B cell multiplets, we filtered out cells that express 408 CD19 in addition to one of the following i) a T cell marker (CD3, CD4, CD8); ii) an innate cell marker (CD14, CD16, 409 CD56). For CD14 + monocyte multiplets, we removed cells that express CD14 and i) a T cell marker (CD3, CD8, CD28); 410 ii) NK marker CD56; iii) a B cell marker (CD19, CD20). For NK cell multiplets, we removed cells that are positive for 411 CD56 and one or more of the following i) T cell markers (CD3, CD4, CD28); ii) monocyte markers (CD14, CD11c); iii) B 412 cell markers (CD19, CD20). For CD4 + T cell multiplets, we removed cells that express CD4 markers (CD3, CD4) and one 413 or more of the following i) CD8 (for CD8 + T), CD56 (for NK), CD19, CD20 (for B cells), CD14, CD16 (for monocytes/NK). 414 Lastly, for CD8 + T cell multiplets, we removed cells that co-expressed CD8 markers (CD3, CD8) and one or more of the 415

Differential expression 431
To identify cell-type-specific response genes and proteins, we conducted pairwise Wilcoxon Rank Sum tests between each 432 PBMC cell-type under stimulation and baseline conditions (e.g., LPS stimulated monocytes vs. baseline monocytes). In 433 each differential expression analysis, genes detected in at least 10% of each cell group were considered. Features with an 434 absolute difference in expression > log2(0.25) and false discovery rate < 5% were regarded as significant results. 435 Benchmarking and inference of single cell pseudotime trajectories were performed with the R package dyno (version 0.1.1) 438 (61). Briefly, for each cell-specific response (e.g., monocytes induced by LPS) trajectories were inferred using Slingshot 439 (the method determined to minimize error and maximize scalability and accuracy). Pseudo-temporal ordering of cells were 440 derived using the corresponding baseline and activated cells as well as the genes and ADTs previously determined to be 441 differentially induced (FDR < 5%) between the two groups of cells. Genes associated with inflammation were obtained from NanoString Technologies, Inc.
(cite, 452 https://www.nanostring.com/products/gene-expression-panels/gene-expression-panels-overview/ncounter-inflammation-panels) 453 and used to calculate an inflammation score for PBMCs at baseline and LPS-stimulated conditions. Feature counts for each 454 cell were divided by the total counts for that cell, multiplied by a scale factor of 10,000, and then log normalized. Normalized 455 read counts for all inflammation-associated genes were divided by the normalized counts of all transcripts and then min-456 max scaled between 0 and 1. Gene lists (Table S3) were used to score cytotoxicity, MAIT or senescence expressions across 457 T and NK cell clusters. To do so, we calculated the mean expression for each cell, within each cluster.

Conflict of Interest 545
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could 546 be construed as a potential conflict of interest.

Acknowledgements 558
We thank all members of the Stitzel and Ucar labs for constructive feedback and discussion on this study and manuscript. 559 We also thank members of the Human Cell Atlas data wrangler team for assistance with depositing the data associated with 560 this manuscript. 561 Table S1: Barcoding information for the 39 antibodies (CITE-seq panel) and cell hashing oligonucleotides (HTOs) used to 563 multiplex treatment information. 564 Table S2: Sequencing quality control information for RNA-seq, cell hashing, and protein quantification data. 565 Table S3: Reference gene lists used to compute single cell scores for inflammation, NK/cytotoxicity, MAIT, TEMRA, and 566 senescence. 567 Table S4: Lists of anti-CD3/CD28 response genes/proteins that are shared or unique to T and NK cell subsets. 568 Table S5: Differential gene and protein expression analysis results for each cell-type at the specified condition vs. that at 569 baseline condition. 570 Table S6: Pathway enrichment analyses for monocyte responses to LPS. 571 Table S7: Genes induced in each of the baseline and LPS-responsive monocyte states (LPS1, LPS2, and LPS3). 572