Single cell multiomic analysis of T cell exhaustion in vitro.

T-cell activation is a key step in the amplification of an immune response. Over the course of an immune response, cells may be chronically stimulated, with some proportion becoming exhausted; an enormous number of molecules are involved in this process. There remain a number of questions about the process, namely: (1) what degree of heterogeneity and plasticity do T-cells exhibit during stimulation? (2) how many unique cell states define chronic stimulation? and (3) what markers discriminate activated from exhausted cells? We addressed these questions by performing single-cell multiomic analysis to simultaneously measure expression of 38 proteins and 399 genes in human T cells expanded in vitro. This approach allowed us to study -with unprecedented depth-how T cells change over the course of chronic stimulation. Comprehensive immunophenotypic and transcriptomic analysis at day 0 enabled a refined characterization of T-cell maturational states and the identification of a donor-specific subset of terminally differentiated T-cells that would have been otherwise overlooked using canonical cell classification schema. As expected, activation downregulated naïve-cell markers and upregulated effector molecules, proliferation regulators, co-inhibitory and co-stimulatory receptors. Our deep kinetic analysis further revealed clusters of proteins and genes identifying unique states of activation, defined by markers temporarily expressed upon 3 days of stimulation (PD-1, CD69, LTA), markers constitutively expressed throughout chronic activation (CD25, GITR, LGALS1), and markers uniquely up-regulated upon 14 days of stimulation (CD39, ENTPD1, TNFDF10); expression of these markers could be associated with the emergence of short-lived cell types. Notably, different ratios of cells expressing activation or exhaustion markers were measured at each time point. These data reveal the high heterogeneity and plasticity of chronically stimulated T cells. Our study demonstrates the power of a single-cell multiomic approach to comprehensively characterize T-cells and to precisely monitor changes in differentiation, activation, and exhaustion signatures during cell stimulation.

INTRODUCTION

The heterogeneity of T‐cells is remarkable; many genes and proteins (“markers”) are associated with cell maturity, trafficking, activation, and function, and these markers can be dramatically modulated over the course of an immune response [1, 2, 3]. Many are dysregulated by the tumor microenvironment [4], so understanding their expression patterns provides critical insight(s) into biological mechanisms of disease, as well as information about potential drug targets. Moreover, it is likely that expression patterns of these markers may be useful in predicting disease, treatment outcome, or therapy‐related adverse events [5].

Technologies to measure T‐cell associated markers have evolved dramatically in the past decade [6, 7]. The use of platforms that assay cells in bulk (like microarrays) has fallen out of favor, because these approaches average expression across many cells, even though the cells vary individually in expression. Single cell RNA sequencing (sc‐RNAseq) overcomes the limitations of bulk measurements, and is powerful because of the large number of transcripts that can be interrogated simultaneously [8]. However, transcription is a “noisy” process that occurs in bursts at irregular intervals and varies even across isogenic/clonal cells [9]; moreover, capture of mRNA in single cell sequencing assays can also be inefficient, leading to drop‐outs of some cell/gene signals. Therefore, cell populations often cannot be clearly resolved based on transcriptional analysis alone.

Cell subset discrimination is greatly enhanced, however, when protein and transcript analysis are combined [10, 11]. Molecular cytometry, a new class of single cell technologies, adapts next generation sequencing (NGS) to single cell analysis to simultaneously provide information about cellular transcripts and proteins. Like flow cytometry, cells are stained with antibodies and unbound antibodies are washed away before analysis. However, unlike flow cytometry, cells are labeled with oligonucleotide‐tagged antibodies (rather than fluorescent molecules), and loaded onto an instrument that captures single cells and lyses them. Beads capture cellular mRNA as well as oligonucleotides associated with cell‐bound antibodies via poly A‐oligo dT interactions. Single cell sequencing then reveals the number of antibodies bound and the target of each bound antibody (as represented by a unique oligonucleotide tag). In addition, cellular transcripts can be examined using targeted panels [12] or whole transcriptome amplification [13]. In this manuscript, we demonstrate the use of molecular cytometry to measure the expression of 38 proteins and 399 targeted T cell transcripts that are relevant to immune cell biology, including T cell exhaustion.

Molecular cytometry technologies, like CITE‐Seq [10], REAP‐Seq [11] and AbSeq [14, 15] carry a number of advantages over other single cell technologies. First, molecular cytometry (also known as genomic cytometry) platforms measure many more parameters than fluorescence or mass cytometry at the single cell level, allowing exquisitely detailed and comprehensive analysis of immune responses. Second, fluorescence and mass cytometry both require subtraction of signals that overlap across channels (compensation) [16, 17, 18]; and this requirement poses significant challenges in terms of experimental design and data analysis. The oligonucleotide tags used in molecular cytometry are unique, which eliminates the need for cross‐channel correction. Third, molecular cytometry can simultaneously interrogate cellular mRNA and proteins more easily than flow or mass cytometry, allowing study of posttranscriptional regulation of protein expression. In sum, molecular cytometry technologies offer deeper profiling of immune responses, as recently demonstrated by two pioneering studies validating the multiomic (BD Rhapsody and BD AbSeq) approach for a comprehensive characterization of either resting or activated immune cells [14, 15].

Upon antigenic challenge, immune responses are maintained and amplified by activated T‐cells, which express unique transcriptional and protein signatures [19]. In the case of an acute infection, markers associated with T cell activation drive key cellular processes, including proliferation, recruitment, homing, cytokine secretion, and cytotoxicity, ultimately resulting in resolution of the immune insult. However, upon sustained antigenic stimulation, as in the case of chronic viral infection or cancer, some activated T cells progressively develop an exhausted phenotype, which is characterized by reduced or lost effector function (e.g., loss of cytokine secretion) and impaired proliferation [20]. During the course of an immune response, the degree to which these exhausted T‐cells are generated dictates key outcomes, including whether a pathogen is cleared or an organism is chronically infected [21], the potential for responsiveness of a tumor to checkpoint inhibition therapy [22], and potentially, the outcome of adoptive immunotherapy using chimeric antigen‐receptor (CAR)‐T cells [23].

The study of activated and exhausted T‐cells offers a unique setting to explore the utility of molecular cytometry‐based immune‐profiling. Further, there remain fundamental gaps in our understanding of what happens when T‐cells are chronically stimulated, namely: (1) what degree of heterogeneity and plasticity do T cells exhibit during chronic stimulation? (2) how many unique cell states (based on transcriptional and protein expression profiles) define chronic stimulation? and (3) what markers discriminate activated from exhausted cells?

We sought to answer these questions by using molecular cytometry to study an in vitro model system of chronic T‐cell stimulation. Our study revealed with unprecedented depth how T cells change upon chronic stimulation at the genotypic and phenotypic level, and also highlighted the use of molecular cytometry for immune profiling, and may identify markers important to study in other settings of chronic T‐cell stimulation, like the tumor microenvironment.

RESULTS

Chronic stimulation of CD8 + T cells in vitro recapitulates phenotypic and functional features of exhaustion

We have developed two in vitro models mimicking either chronic or transient T‐cell stimulation. Chronic stimulation was achieved by continuously stimulating T cells with recombinant human IL‐2 (rhIL‐2) and αCD3/CD28 beads for 14 days. Transient stimulation was achieved by stimulating T cells with rh‐IL2 and αCD3/CD28 beads for 3 days, followed by resting in the presence of rhIL‐2 for additional 11 days. Cells were collected at different time points, as indicated in Figure 1A.

FIGURE 1

In vitro chronic stimulation recapitulates features of T‐cell exhaustion. (A) Depiction of the in vitro model used for chronic and transient stimulation of T cells isolated from healthy donors. Cells were collected and frozen at day 0, 3, 7, and 14 for downstream analysis. (B) Representative flow cytometry analysis of the expression of inhibitory receptors CD279 (PD‐1) and CD223 (LAG‐3), as well as cytokines (IFNg, TNF, and IL2) on CD8+ T cells in response to chronic or transient stimulation. (C) Representative flow cytometry analysis of CD4+ T‐cells. Analyses performed on fresh cells from at least three independent experiments

To assess whether CD8+ T cells in our in vitro model system acquired phenotypic and functional features characteristic of chronically stimulated T cells, we measured the upregulation of well‐characterized inhibitory receptors, and the production of inflammatory cytokines using flow cytometry. CD8+ T cells demonstrated appreciable upregulation of the inhibitory receptors CD279 (PD‐1) and CD223 (LAG‐3) upon 3 days of stimulation (Figure 1B). LAG‐3 expression was maintained through day 14 upon chronic stimulation, with a gradual downregulation of PD‐1 for some cells, as previously described. Under transient stimulation conditions, CD8+ T cells progressively lost expression of both inhibitory receptors by day 14. Next, we analyzed T‐cell function at each time point, and for each experimental condition (i.e., chronic vs. transient stimulation). At baseline and after 3 days of culture with rhIL‐2/αCD3/CD28, cells produced IFNγ and/or IL2; however, at day 7–in both chronic and transient conditions—cytokine expression was greatly reduced (Figure 1B). This observation suggests that cells become functionally impaired beyond 3 days of stimulation. Cytokine production after 14 days was dramatically impaired when cells were stimulated continually (chronic stimulation), while T‐cell function was recovered in transient stimulation conditions (Figure 1B). Similar results were observed for CD4+ T‐cells (Figure 1C). Together, these data demonstrate that well‐established phenotypic and functional changes associated with T‐cell exhaustion [20] are specifically and gradually induced in our in vitro model system.

AbSeq enables protein detection with specificity and resolution comparable to flow cytometry

Having confirmed that our in vitro system models the dynamics of marker expression associated with T cell exhaustion, we next compared the expression of different cell surface proteins, as measured by flow cytometry and AbSeq. The flow cytometry and AbSeq panels used for the comparison are outlined in Table 1. Figure 2A shows CD39 expression on total T cells (CD8+ and CD8‐ [mostly CD4+]) over the time course of the study. CD39 expression on CD8+ and CD8‐ T‐cells was detectable in a small subset of cells upon 7 days of chronic stimulation. After 14 days of stimulation, the frequency of CD8+ and CD8‐ T cells expressing CD39 increased, along with the level of CD39 expression. In contrast, no CD39 expression was detected in any T‐cell subset after 14 days of transient stimulation. Once detectable, CD39 expression patterns were qualitatively similar for both AbSeq‐ and flow cytometry‐based measurements. For 7 of 10 proteins whose detection was compared using the two approaches, AbSeq and flow cytometry reported similar frequency of cells expressing each marker (Figures 2B and S1A). Additionally, the kinetics of marker expression were relatively consistent between the flow cytometry and AbSeq measurements (Figures 2C and S1); however, sample size for this analysis was limited. Our data suggest that, consistent with past reports, sensitivity and specificity of the AbSeq approach are generally comparable to flow cytometry.

TABLE 1

Multiplex panels used for protein analysis. A 13‐parameter panel was used for flow cytometry‐based analysis. A 38‐parameter panel was used for AbSeq‐based analysis. Twelve specificities were common to both panels and were used to compare the performance of the two platforms

Marker	Clone	Flow cytometry or AbSeq	Fluorochrome/dye used for flow cytometry analysis
CD4	SK3	Both	BUV805
CD8	RPA‐T8	Both	BUV395
CD45RA	HI100	Both	APC‐H7
CD62L	DREG‐56	Both	FITC
CD95	DX2	Both	BV786
CD279 (PD‐1)	EH12.1	Both	PE‐Cy7
CD223 (LAG‐3)	T47‐530	Both	BV480
CD366 (TIM‐3)	7D3	Both	BV711
CD357 (GITR)	V27‐580	Both	BV421
CD152 (CTLA‐4)	BNI3	Both	PE
CD39	TU66	Both	BUV737
CD103	Ber‐ACT8	Both	APC
Live/Dead	N/A	Flow Cytometry Only	7‐AAD
CD3	SK7	AbSeq Only	N/A
CD14	MфP‐9	AbSeq Only	N/A
B7‐H4	MIH43	AbSeq Only	N/A
CD127	HIL‐7R‐M21	AbSeq Only	N/A
CD134	ACT35	AbSeq Only	N/A
CD137	4B4‐1	AbSeq Only	N/A
CD154	TRAP1	AbSeq Only	N/A
CD183	1C6/CXCR3	AbSeq Only	N/A
CD185	RF8B2	AbSeq Only	N/A
CD194	1G1	AbSeq Only	N/A
CD196	11A9	AbSeq Only	N/A
CD197	3D12	AbSeq Only	N/A
CD2	RPA‐2.10	AbSeq Only	N/A
CD25	2A3	AbSeq Only	N/A
CD27	M‐T271	AbSeq Only	N/A
CD270	CW10	AbSeq Only	N/A
CD278	DX29	AbSeq Only	N/A
CD28	CD28.2	AbSeq Only	N/A
CD30	BerH8	AbSeq Only	N/A
CD38	HIT2	AbSeq Only	N/A
CD49a	SR84	AbSeq Only	N/A
CD54	HA58	AbSeq Only	N/A
CD69	FN50	AbSeq Only	N/A
CD7	M‐T701	AbSeq Only	N/A
CD94	HP‐3D9	AbSeq Only	N/A
CD98	UM7F8	AbSeq Only	N/A

Note: “N/A” = not used.

FIGURE 2

AbSeq and flow cytometry enable protein measurement with equivalent specificity and sensitivity. (A) Qualitative analysis of CD39 kinetic expression within CD8+ and CD8− subsets of CD3+ T cells performed using either flow cytometry (top panel) or AbSeq (bottom panel). (B) Percentage of CD8+ cells expressing PD‐1, GITR, TIM‐3 and CD39 at day 0 (D0), day 3 (D3), day 7 (D7 C) and day 14 (D14 C) of chronic stimulation using either flow cytometry (red line) or AbSeq (blue line). (C) Levels of expression of PD‐1, GITR, TIM‐3 and CD39 on CD8+ T cells at day 0 (D0), day 3 (D3), day 7 (D7 C) and day 14 (D14 C) of chronic stimulation. Mean fluorescence intensity (red line, left y‐axis) and mean molecular count (blue line, right y‐axis) were used to measure relative antigen expression levels using flow cytometry and AbSeq, respectively. This side‐by‐side analysis was performed on different aliquots of the same cells derived from the same donor (Donor 1), cryopreserved at the indicated time points. Flow cytometry data were down‐sampled in order to analyze the same number of cells across the two platforms at every time point (range 528–2446 cells)

Molecular cytometry enables deeper profiling of effector T cells

Next, we conducted a multi‐omic analysis of resting CD8+ T cells (day 0) by simultaneously measuring the expression of 38 proteins and 399 T cell‐specific genes at the single cell level. The 38‐plex AbSeq panel is outlined in Table 1; and the list of the 399 genes measured by the targeted RNAseq panel is reported in File S1. First, we used AbSeq to quantify the different subsets of unstimulated CD8+ T cells (day 0) based on differential expression of markers commonly used for identification of T‐cell differentiation states [24, 25]. CD45RA and CD28 were used to define CD45RA− CD28+ central memory cells (CM), CD45RA− CD28− effector memory cells (EM) and CD45RA+ CD28− effector memory RA cells (EMRA). CD27 was additionally used for the identification of CD45RA+ CD28+ CD27high naïve cells (N). As expected, both qualitative and quantitative analysis revealed differences across the three donors when looking at the different T cell subsets (Figure 3A,B). We also noted a unique population of CD45RA+CD28+CD27low cells in Donor 3, which represented ~12% of the CD8+ T cell compartment.

FIGURE 3

Deep characterization of fresh CD8+ T‐cell maturational states. (A) Gating strategy used to identify CD8+ naïve (blue box), central memory (CM, purple box), effector memory (EM, yellow box), effector memory RA (EMRA, red box) T cells from three donors at day 0 based on measurement of CD45RA, CD28 and CD27 expression via AbSeq. A unique population of CD45RA+CD28+CD27low cells (CD27low, green box) was detected in donor 3. (B) Frequency of CD8+ T‐cell subsets across the three donors. (C) Single‐cell heatmap of selected proteins (Ab) and genes differentially expressed across the three donors (fold change ≥2, q ≤ 0.05). Four hundred cells per donor are represented. Each column represents a single cell from total CD8+ cell gate. Event columns are colored on a 0–100% min‐max pseudocolor scale based on relative parameter expression. Markers are grouped according to the cell type they are commonly expressed on, with CD45RA repeated because it is expressed by both naïve and effector memory cells; for Donor 3, the heatmap values reflect mostly contributions from effector memory cells. (D) t‐SNE visualization of the CD8+ T‐cell subsets showing different cell clusters with a naïve, CM, EM, EMRA and CD27low/− phenotype across the three donors. T‐SNE plots were generated based on expression of 38 proteins and highly dispersed genes. Cells were color‐coded based on the phenotype described in panel A. The red dashed circle indicates EMRA and CD27low/− cells clustering together in Donor 3

To further define proteins and transcripts associated with the different maturation states across donors, we performed differential gene and protein expression profiling comparing naïve cells to CM, EM, and EMRA T cells. Consistent with published findings [25], we observed higher levels of CD45RA, CD62L, CD197 (CCR7) proteins and SELL, LEF1, IL7R and CCR7 transcripts in naïve cells (File S2). Also, as expected, the multi‐omic analysis revealed that EM and EMRA cells (yellow and red boxes) express high levels of PD‐1 protein and transcripts encoding cytotoxic proteins NKG7 (NKG7), granzymes (GZMA, GZMB, GZMH), granulysin (GNLY), and perforin (PRF1) (File S2). Interestingly, CD8+ T cells from Donor 3 showed higher expression of effector‐associated markers, compared to samples from Donors 1 and 2 that were enriched for cells displaying a naïve/central memory profile (Figure 3C). (Note that because CD45RA is a marker expressed on both naïve and effector cells, the data for this marker are repeated in the heat map.) Also, our analysis revealed new markers not previously associated with the different T cell maturation states (PIK3IP1, PASK, and TXK in CD8+ N cells and DUSP1 and IFNGR1 in CD8+ CM cells).

To better characterize the CD45RA+CD28+CD27low population in Donor 3, we generated t‐distributed stochastic neighbor embedding (t‐SNE) maps for each donor. The events in the t‐SNE plot are color‐coded based on the maturation state of cells as defined by differential expression of CD45RA, CD28 and CD27 (Figure 3D). The high number of EMRA cells in Donor 3 formed a distinct cell cluster that included the CD27low cells (Figure 3C, bottom panel, red dashed circle), suggesting that CD27low cells share gene and protein expression patterns with EMRA cells rather than CD45RA+CD28+CD27high naïve T‐cells, as might have been expected. To test whether CD27low cells are more closely associated with EMRA than naïve cells, we performed differential expression analysis, comparing all measured transcripts and proteins between CD27low and naïve cells derived from Donor 3. We found that genes and their corresponding proteins common to naïve cells (and absent from EMRA), like SELL/CD62L, IL7R/CD127, CCR7/CD197, were expressed at higher levels in naïve cells compared to CD27low cells (Table S1). Conversely, CD27low cells were enriched for proteins and genes expressed by EMRA cells, like PD‐1, PRF1, GZMB, and IFNG (Table S2). Notably, 19 of the 20 markers detected at higher levels in CD27low cells (compared to naïve cells) were also enriched in EMRA cells (Table S3). In sum, our results demonstrate that CD27low cells, that might have been wrongly classified as naïve cells based solely on expression of CD45RA, CD28, and CD27, are in fact, likely to be a subset of effector memory cells capable of reexpression of both CD45RA and CD27. This analysis reveals the power of multi‐omic (protein and mRNA) analyses for more precise identification of cell types, and for deeper profiling of classical T‐cell phenotypes.

Molecular cytometry identifies gene and protein signatures associated with different modes of T cell activation

Having confirmed the utility of molecular cytometry for detailed cellular profiling, next, we used this approach to investigate temporal changes in protein and gene expression that accompany chronic T‐cell stimulation and exhaustion. The same AbSeq and targeted gene panels (Table 1 and File S1) that were used for the characterization of T‐cell maturational states, were used for a comprehensive analysis of CD8+ T cells at different stages of chronic and transient stimulation. High‐dimensional datasets were generated for each donor at the different time points based on expression of 38 proteins and highly dispersed genes. Projection of cells from the of the three concatenated donors using t‐SNE revealed three major cell groups (Figure 4A, “all donors”). The first group encompassed events from resting cells (before activation, “day 0,” red colored events), while the second group included events from 3, 7, and 14 days of chronic stimulation (blue, orange, and green colored events), and the third group comprised cells from the 3‐day stimulation +11‐day rest condition (“Day 14 Transient,” purple colored events). These results demonstrate that expression changes among the 38 markers measured by AbSeq in conjunction with targeted gene expression analysis clearly capture the distinction between resting and activated cells. Notably, the combined analysis of mRNA and protein provided the best discrimination of cell subsets at stimulation time points (Figure S2). For example, when transcripts alone are used to discriminate stimulation time points, a population of cells stimulated and then rested‐down (day 14 Transient, purple, mRNA only; see black arrow) appears to track with unstimulated cells (red); however, when mRNA and protein are combined, no cells have overlapping profiles at these time points (Figure S2; mRNA+Ab). Similarly, when antibodies alone are examined (Ab Only) a population of cells stimulated for 3 days falls in a contiguous group with cells stimulated for 7 days (black arrow). However, when mRNA and protein analysis are combined, these cells are mostly distinct. Thus, a multiomic approach appears to provide better discrimination of resting and activated CD8+ T‐cell subsets.

FIGURE 4

Identification of common signatures associated with T‐cell activation. (A) t‐SNE visualization of CD8+ cells clusters at day 0 (red), day 3 (blue), day 7 and 14 of chronic stimulation (orange and green, respectively), and day 14 of transient stimulation (purple). t‐SNE plots were generated based on expression of 38 protein and highly dispersed genes. (B) Venn diagrams indicate the number of shared or uniquely upregulated genes and proteins across the three donors at day 3, day 7 chronic stimulation (D7 C), day 14 chronic stimulation (D14 C), and day 14 transient (D14 T) stimulation, as compared to unstimulated cells (day 0)

We also examined CD4+ T‐cells in our model system, and found that a number of transcripts and proteins were elevated in both CD4+ and CD8+ T‐cells chronically stimulated for 14 days (compared to resting cells; File S3). These included a chemokine (CCL3), cytokines (IFNg, CSF2), a marker of cellular activation (CD25), cytotoxic enzymes (granulysin, granzyme B), markers of exhaustion (CD357, LAG3), transcription factors (TYMS, UBEC2C, and ZBED2), and a marker of cell proliferation (PCNA). These molecules provide a signature, common to CD4+ and CD8+ T‐cells, for activated cells.

To identify the genes and proteins uniquely associated with the mode of T‐cell activation (chronic versus transient), we performed differential expression analysis for each individual donor by comparing each time point of the in vitro activation system to each other (File S4). By comparing the individual lists of differentially expressed genes and proteins, we identified, for example, common signatures defined by markers upregulated in each donor at each time point, as compared to unstimulated cells (day 0). The results of this analysis are summarized using the Venn diagrams in Figure 4B. We also observed upregulation of genes and/or proteins shared by two out of three donors, or unique to each donor. The complete list of shared and unique markers upregulated at each time point (in each donor), as compared to day 0, is reported in File S5.

Next, we focused on genes and proteins that were upregulated 3‐fold or higher upon activation in at least two of the three donors. This approach resulted in the identification of four sets of genes and proteins with unique expression patterns. The first set (Figure 5A, red bar) consisted of markers elevated after 3 days of stimulation (D3), but downregulated thereafter (CD278, CD69, IFNg, IL9, and Lymphotoxin A [LTA]). The second set (Figure 5A, blue) represented markers elevated after 14 days of chronic stimulation (D14C), and included TNFSF10, YBX3, CSF2, BIRC3, ENTPD1, and CD39. The third set of markers (Figure 5A, green) was upregulated at all stimulation time points, but generally downregulated when the stimulus was removed (as measured after 3 days stimulation and 11 days rest; D14T in Figure 5A). This set included genes and proteins that could be broadly classified into markers of activation/proliferation (CD25, CD357, CD54, CD98, CD137, GZMB, IL2Ra, PCNA, TOP2A, and TYMS) versus inhibition/exhaustion (CD223, IRF4, LAG3, LGALS1, and ZBED2). The fourth set (Figure 5A, purple) encompassed markers whose expression was upregulated in cells transiently stimulated. Taken together, our analysis reveals unique activation signatures that can be differentially linked to the duration of activation.

FIGURE 5

Identification signatures associated with distinct modes of T‐cell activation. (A) Markers whose expression is upregulated (≥ three‐fold, in at least two donors; left y‐axis), as compared to day 0, at each time point (right y‐axis). Markers uniquely upregulated with 3 days of stimulation are indicated by the red bar, while markers unique to chronic 14‐day stimulation are indicated by the blue bar. The green bar indicates which markers are elevated at all time points in the chronically stimulated cells, while the purple bar denotes markers uniquely associated with the 3‐day stimulation/11 day rest (transient stimulation, D14T) condition. (B) Single‐cell heatmap of proteins (Ab) and genes upregulated in all three donors at day 0 and downregulated upon cell activation (fold change ≥2, q ≤ 0.05). Data from 420 cells measured at each time point. An equal number of cells (140) from each of the three donors is represented for each time point

Notably, our analysis also revealed a set of genes and proteins that were selectively associated with resting cells in all three donors. These markers were down‐regulated upon stimulation, and then partially reacquired after the stimulation was removed. This marker set included proteins and genes associated with naïve T‐cells (CD45RA, IL7R, and its corresponding protein CD127), along with DUSP1, FOSB, and JUN (Figure 5B).

Molecular cytometry reveals unique relationships between inhibitory and proliferation markers

To understand whether markers of activation or inhibition are expressed on distinct cells or if these genes and proteins could be co‐expressed on the same cells, we plotted expression in bivariate plots using concatenated data from the three donors. Among the various combinations of markers, those involving LGALS1 and ZBED2 transcripts stood out for their progressive expression over the time course. At baseline, very few cells expressed LGALS1 (the mRNA for Galectin‐1, an immune inhibitory molecule) or proliferation markers TYMS and PCNA (Figure 6A). At day 3, LGALS1 and both proliferation markers were upregulated, and cells expressing all possible combinations of markers (i.e., one marker alone, neither marker, both markers) could be detected (Figure 6A). However, after 7 days of stimulation, almost all cells expressing TYMS or PCNA co‐expressed the inhibitory molecule LGALS1, suggesting the inhibitory potential of these cells. Conversely, after only 3 days of stimulation, the great majority of cells expressing TYMS and PCNA co‐expressed ZBED2, a transcription factor associated with progression to T‐cell dysfunction (Figure 6B). Beyond 3 days, the frequency of cells expressing ZBED2 was gradually reduced, with a concurrent increase in the frequency of cells single positive for PCNA and TYMS (Figure 6B). Thus, expression of LGALS1 and ZBED2 appear reciprocal between days 3–14 (Figure 6C) of stimulation.

FIGURE 6

Combinatorial analysis of inhibitory and proliferation marker expression. Bi‐variate plots showing the relationship between the inhibitory markers LGALS1 (A) or ZBED2 (B) and the proliferation markers TYMS or PCNA throughout chronic stimulation. (C) Frequency of CD8+ T‐cells expressing LGALS1 and ZBED2 throughout chronic stimulation. Data were generated from concatenated samples from three donors

Comparing chronic to transient stimulation

We also compared transient stimulation (3 days of stimulation, followed by 11 days of rest) to chronic stimulation (14 days of continual stimulation). Differential gene/protein analysis revealed that chronically stimulated cells are enriched for a particular set of markers (compared to transiently stimulated cells, Table 2), which includes markers of cytotoxic cells (GZMB, GNLY, CD94, GZMH, TNFSF10), suppression (CD39, LAG3, CD223), effector function (CSF2 (GMCSF), CCL3 (MIP1α), IFNg, CCL4 (MIP1β)), and immune checkpoints (CD357, CD137). In contrast, transiently stimulated cells appeared capable of retaining (or resting back down to) a naïve‐like phenotype after the initial 3‐day stimulation, as evidenced by higher levels of CD45RA, CD28, CD27, CD127, and IL7R (Table 2). The single cell data show that fold‐change (a bulk measurement) is accurately reflecting the single cell co‐expression of these molecules; thus, the cells we are detecting at day 14 in the transient condition are naïve cells, not central memory cells (which can share expression of some of the molecules reported above, data not shown). An alternate interpretation of this data is not that cells rest back down to a naïve‐state, but rather that transiently stimulated cells died, leaving mostly naïve cells that were capable of persisting in the culture conditions. In any case, cells transiently stimulated and then rested are distinct from unstimulated (data not shown) and from cells after 3 days of stimulation (Table 3).

TABLE 2

Genes and proteins (ab) that are differentially expressed between cells chronically stimulated for 14 days versus those stimulated for 3 days and then rested for 11

Elevated in D14 chronic versus transient		Elevated in D14 transient versus chronic
Marker	Fold change	Marker	Fold change
GZMB	21.0	CD45RA (Ab)	8.3
CSF2	13.5	CD28 (Ab)	5.0
GNLY	13.4	CD27 (Ab)	4.7
CD94 (Ab)	11.1	IL7R	4.1
CD357 (Ab)	10.9	CD183 (Ab)	3.8
CCL3	10.7	CD44	3.3
IFNG	10.6	TRIB2	3.0
CCL4	10.3	S100A10	3.0
ZBED2	9.4	CD3 (Ab)	2.9
CD39 (Ab)	7.1	CD127 (Ab)	2.5
CD25 (Ab)	6.8	FAM65B	2.4
CD54 (Ab)	6.8	MYC	2.4
LAG3	5.8	PIK3IP1	2.4
CD223 (Ab)	5.5	BIN2	2.3
DUSP4	4.4	GZMA	2.2
CCL1	4.4	CD8(Ab)	2.2
GZMH	4.1	CCR2	2.1
TNFSF10	3.9	CCR5	2.1
KLRC1	3.7	CD95 (Ab)	2.1
CD137 (Ab)	3.6	CCL5	2.1

TABLE 3

Genes and proteins (ab) that are differentially expressed between cells transiently stimulated for 3 days then rested for 11 versus cells stimulated for 3 days

Elevated in D14 transient over D3
Marker	Fold change
CD52	6.2091
CD62L (Ab)	5.0631
CD27 (Ab)	5.0133
CD38 (Ab)	4.3823
CD183 (Ab)	4.0164
SELL	3.7154
KLRK1	3.4889
TRIB2	3.3777
IL7R	3.2278
CD8 (Ab)	3.0546
GZMA	3.0258
SELPLG	2.9951
IFITM2	2.7836
PIK3IP1	2.7315
BIN2	2.6601
FYB	2.6171
FAM65B	2.3680
IFITM3	2.3429
XBP1	2.2892
CD127 (Ab)	2.2744
CD27	2.2338
ARL4C	2.2228
CD44	2.1958
DPP4	2.1885
CD37	2.1157
CCR2	2.1125
ITGA4	2.1097
S100A10	2.0985
GIMAP2	2.0716
LEF1	2.0246
ZNF683	2.0114
LGALS3	2.0110
TSPAN32	2.0048
ADA	2.0041

We next used Phenograph to define clusters of cell types/states present within the transiently‐stimulated condition. Although these cells mostly‐expressed markers of naïve cells, there were clusters of cells that retained signatures of activated/proliferating cells. One cluster expressed CCL3, CCL4, IFNg, and CD69 for example, while another expressed TOP2a, PCNA, and MKI67 (data not shown). These results reveal the heterogeneity associated with length of stimulation, and illustrate how the combination of AbSeq and Phenograph are powerful tools for cellular profiling.

Interestingly, when we analyzed TCF7 expression through the stimulation time course, we found a difference between TCF7+ and TCF7‐ cells arise after 14 days of stimulation. At this time point, we found that, compared to PD1+ TCF7‐ cells, PD1+ TCF7+ cells were enriched for: CCR7 and CD27 (characteristic of naïve and central memory T‐cells; 2.2 and 2.0‐fold enrichment respectively), CXCR3 (a marker of naïve CD8+ T‐cells capable of enhanced differentiation into effector cells [reference]; 1.6‐fold enrichment), and CXCL13 (which is thought to enhance development of tertiary lymphoid structures, and may help convert Treg‐mediated immune suppression to de novo activation of an adaptive immune response; 7.1‐fold enrichment). Thus, PD1+ TCF7+ cells preferentially express markers that are associated with the potential for maintenance and amplification of immune responses; these findings support recent reports that PD1+ TCF7+ expression may denote cells that can be rescued from immunosuppression.

Molecular cytometry reveals correlation between Cellular Transcription & Translation

In addition to detailed cellular profiling, molecular cytometry technologies provide a unique ability to correlate the expression of genes and their corresponding proteins at the single‐cell level. In this study, we examined this correlation for 23 pairs of genes and their corresponding proteins, and determined if kinetic changes (over the activation time course) were similar. In some cases, mRNA and the corresponding protein exhibited concordance in terms kinetics of expression. For example, ENTPD1 and its corresponding protein CD39 were concordantly upregulated at day 14 of chronic stimulation (Figure 7A). Similarly, both IL7R and its corresponding protein CD127 were concordantly downregulated upon chronic cell stimulation. For certain markers, we observed discordant patterns for mRNA and protein expression. For example, unstimulated cells at day 0 express basal levels of CD69 mRNA but not protein. Upon 3 days of stimulation, we observed CD69 mRNA downregulation, with concomitant upregulation of CD69 protein. After 3 days of stimulation, both CD69 mRNA and protein were downregulated. Finally, in some cases, protein expression was detected (and varied over the course of the stimulation), with rare and low mRNA expression (PDCD1 mRNA/PD1 protein, Figure 7A; other cases in Figure S3). In sum, the relationship between protein and mRNA expression was complex, and observed in all possible patterns (concordant/discordant in frequency of expression; concordant/discordant in level of expression). This is consistent with previous reports describing the stochastic nature of mRNA transcription [26].

FIGURE 7

Correlation between gene and protein expression at the single cell level. (A) Kinetic analysis of mRNA and protein expression. mRNA and protein levels were measured as mean molecular count (MMC) and depicted by the red trace (protein) and blue trace (mRNA) at day 0 (D0), day 3 (D3), day 7 (D7 C) and day 14 (D14 C) of chronic stimulation. Analysis performed on three individual donors. Data represented as mean ± SD. (B) Single‐cell heatmap of gene and protein expression. A single cell is represented in each column. One hundred cells from a single donor (Donor 2) are represented at each time point. Event columns are colored on a 0–100% min‐max pseudocolor scale based on relative parameter expression

The above analyses summarize the data on protein/gene expression correlation, but do not illustrate the heterogeneity of mRNA/protein concordance on a cell‐by‐cell basis. Concordance and discordance at the single cell are shown by heat map in Figure 7B. These heat maps reveal the complexity of mRNA/protein expression and kinetics at the single‐cell level. In accordance with the similar kinetics observed for ENTPD1 and CD39 (both upregulated after 14 days of chronic stimulation; Figure 7A), the great majority of cells (represented by colored bars) co‐expressed mRNA and protein (Figure 7B). The same was not true for IL7R and CD127. Despite the overall similarity in kinetics between IL7R and CD127 (expressed at day 0, reduced at day 3, almost off on days 7 and 14; Figure 7A), a higher degree of heterogeneity was observed, with variable numbers of cells expressing protein only, gene only or both over time. In contrast, CD69 mRNA is expressed by many cells at all‐time points (Figure 7B), but protein is expressed at essentially a single time point (day 3).

DISCUSSION

Molecular cytometry is a potentially powerful method for immune monitoring; however, the technology is relatively new. We sought to provide data that qualifies AbSeq, a new molecular cytometry technology, for use in immune monitoring and to demonstrate the power of the approach. We first provide data suggesting that AbSeq compares favorably with the gold standard flow cytometry. As past studies [10, 11, 14] have shown, the proportion of cells expressing various markers are similar across the two platforms. In fact, when protein expression is visualized using bivariate plots, and as expected when using the same antibody clone for the two different methods, the patterns are quite similar across the technologies. Where we found discordance in the percentage of cells expressing markers such as LAG‐3, CD62L, and CTLA‐4, we speculate that this was a function of poor performance of the fluorochrome‐conjugated antibody, since it is well‐recognized that performance of an antibody clone may depend on the fluorochrome chosen or the design of a multicolor panel [16]. To definitively demonstrate this, a variety of oligo and fluorescence tags would need to be compared for the same antibody clone; nevertheless, the issue highlights two important advantages of molecular cytometry over flow cytometry. First, all antibody tags are oligonucleotide sequences, so they share similar chemical and biophysical properties. As such, antibody performance in molecular cytometry is not dependent on the tag chosen or labeling method. Second, signal, background, and sensitivity are similar for all antibody tags that are sequenced, unlike flow cytometry, where the detectors (i.e., photomultiplier tubes) used for detection of each fluorochrome can vary a great deal [27]. Moreover, in flow cytometry, the signals observed for a given level of protein expression may differ by fluorochrome, whereas the digital nature of molecular cytometry signals allow a relatively consistent relationship to protein expression. In principle, quantification of receptor levels by molecular cytometry may therefore be considerably more straightforward than flow cytometry, although the efficiency with which oligonucleotide‐tagged antibodies are captured onto beads and the depth of sequencing may complicate this application.

Unlike other molecular cytometry publications, we measured kinetics of marker expression, which provide further evidence of the concordance between flow and molecular cytometry. Our kinetic analysis also allowed comparison between protein and transcript expression for the 23 mRNAs with a coordinate antibody target. We found great variety in the relationships between transcript and protein expression among individual markers; we did not find high correlation between mRNA and protein expression, and often observed transcripts that were up‐regulated or present without concomitant protein expression, and vice versa. In some cases, discordance may reflect protein or transcript expression that is below the level of detection. However, we often found—for phenotypically similar cells–that one cell had high levels of transcript while another cell did not. The relationship between transcript and protein expression did not follow any clear rules; it was not specific to proteins or cell types that shared a function, nor were there common kinetic patterns. This suggests that discordance between transcript and protein expression may reflect (a) the stochastic nature of transcript expression, which occurs in bursts [28, 29, 30] or (b) regulation of transcript expression through mechanisms that are specific to individual proteins. In either case, it is very clear that transcript expression is not a direct surrogate for protein expression, nor is the abundance of an mRNA correlated with protein [31]; these factors could diminish the value of platforms that solely measure single cell RNA expression in isolation (without consideration of protein expression) for immune monitoring.

Our analysis of naïve, memory, and effector T‐cell subsets demonstrates the power of molecular cytometry for accurate identification of cell populations. In our oldest study subject (Donor 3, 61‐years‐old), we observed a population of CD45RA+ cells with intermediate expression of CD27. Using canonical cell classification schema [25], such cells might have been categorized as naïve; however, by more deeply profiling cells–measuring the many proteins and genes available in our experiment–it became clear that CD27low cells were much more similar to effector T‐cells than naïve T‐cells. The broader lesson revealed by this example is that molecular cytometry, and other high‐parameter single‐cell technologies can check and challenge the classical schema used to categorize the maturational status of cells. The problem presented by unusual patients or cell populations, which may be outliers in immune monitoring studies, can also be mitigated by molecular cytometry. The highly multiplexed, and simultaneous measurement of protein and mRNA, provides increased context that better describes unusual cell populations and allows more accurate classification and enumeration. This advantage, however, may be mitigated somewhat by the limited throughput, lower cell numbers, and expense of molecular cytometry platforms.

Our study examined, with great depth, how T cells change with activation, using molecular cytometry to analyze time course specimens from an in vitro activation model. The model recapitulated the process of T‐cell activation and exhaustion [24], as confirmed by flow cytometric analysis of checkpoint molecules and cytokines. We identified sets of proteins and genes that are uniquely upregulated (compared to resting cells) at each time point we analyzed. Five markers were upregulated only at the 3‐day time point, including the genes coding for known activation marker CD69 and the effector cytokine interferon gamma. We posit that these markers, whose expression during short‐term, ex vivo stimulation (6–24 h) is well‐documented by flow cytometry [32], remain elevated even after 3 days. Notably, expression of IL9 and lymphotoxin A (LTA) mRNA were also highly upregulated at this time point. These represent new targets for immune assessment, perhaps using intracellular cytokine staining (ICS) by flow cytometry (provided that these mRNA are translated into protein). The identification of new immune monitoring targets–beyond the common cytokines measured by ICS–demonstrates the value of molecular cytometry.

We also identified the proteins and genes upregulated throughout the time course, which fell neatly into two groups—those that are associated with enhancing cell function and those that inhibit cell activity. Unlike bulk assays, in which thousands of cells are averaged for analysis, single cell molecular cytometry data allowed us to ask whether the upregulation of these markers was associated with two distinct cell types (e.g., activated vs. exhausted) or whether these markers could be co‐expressed. We found that many activation and inhibition markers were co‐expressed, suggesting great plasticity in cell state and function, and that the relationships between some of these markers changed over the stimulation period. For example, markers associated with cell proliferation PCNA and TYMS were expressed by cells after 3 days of stimulation, with and without LGALS1 expression. However, by day 7 of our stimulation assays, all cells expressing proliferation antigens co‐expressed LGALS1. This gene transcript is notable because it encodes the Galectin‐1 protein, which is known to inhibit cell proliferation in the tumor microenvironment and is a target of immunotherapy agents [33]. Our result suggests the possible existence of an autocrine or paracrine feedback loop regulating cell proliferation, involving PCNA /TYMS and Galectin‐1. Disruption of this feedback loop, using antibodies to Galectin‐1, may provide a new means to prevent exhaustion of CAR‐T cells during their manufacture. Similarly, the ability of CAR‐T cells to persist in vivo might be enhanced by genetically engineering the over‐expression of markers that are normally downregulated with stimulation, including DUSP1, FOSB, and JUN (which were part of our resting cell signature). Indeed, a recent report describing “exhaustion‐resistant” CAR‐T cells with over‐expressed c‐Jun supports this possibility [34]. Finally, measurements of JUN, LGALS1, or the suite of markers upregulated with 14 days of stimulation may provide a predictor or indicator of the capacity of a CAR‐T cell product to persist in vivo, raising the intriguing possibility of a companion diagnostic for CAR‐T cell therapy.

The high parameter data provided by single cell molecular cytometry offers an unparalleled tool to better define a molecule's function, expression, and disease relevance. For example, our study also reports reciprocal expression of LGALS1 and ZBED2, in relationship to the expression of proliferation markers PCNA and TYMS; expression of LGALS1 is gained in proliferating cells over stimulation, while ZBED2 is lost from proliferating cells over the course of stimulation. ZBED2 expression has recently been shown to mark a subset of melanoma‐infiltrating CD8+ T cells poised to progress to a dysfunctional, exhausted state [35]. The same study also described highest and lowest proliferative potential at early and late stages of exhaustion, respectively, thus corroborating the importance of simultaneously assessing inhibitory and proliferation markers at the single‐cell level.

We have also identified a set of markers exclusively expressed upon 14 days of chronic stimulation. Among these markers are the gene ENTPD1 and its corresponding protein CD39, recently described as a marker defining tumor antigen‐specific, exhausted TIL [36]. The significant loss of cytokine production observed by flow cytometry with chronic stimulation suggests that expression of this set of markers correlates with T‐cell dysfunction. Further studies investigating the expression of the identified signatures of primary tumor infiltrating lymphocytes are required to validate these hypotheses.

Ultimately, the power of molecular cytometry lies in the rich datasets it provides. These can be produced with fewer antibody panel design complications than other cytometry technologies. The technology approach offers an important advance in our ability to characterize–and exploit–cellular immunity.

MATERIALS AND METHODS

Cell preparation and cryopreservation

Peripheral blood was collected from N = 4 healthy donors (age 28‐, 32‐, 36‐, and 61‐years‐old) in accordance with approved IRB protocol BDX‐ASCP0. Peripheral blood mononuclear cells (PBMCs) were isolated using Ficoll‐Paque Plus (GE Healthcare). T cells were isolated using BD IMag™ Human T Lymphocyte Enrichment Set, as per manufacturer's instructions (BD Biosciences). Aliquots of cells collected at different time points of stimulation (day 0, 3, 7, and 14) were cryopreserved in freezing medium containing 90% fetal bovine serum (FBS; Hyclone Laboratories Inc.) and 10% DMSO (Sigma‐Aldrich) and stored in liquid nitrogen. Cryopreserved cells were thawed in a 37°C water bath, diluted with 1 ml of warm complete culture medium, then transferred to a tube containing 10 ml of warm complete culture medium. Cells were centrifuged at 300g for 5 min prior to downstream processing. Cell viability of fresh or thawed cells was overall consistently ≥90%.

T‐cell stimulation

Freshly isolated T cells were plated in a 24‐multiwell plate (Corning) at the concentration of 1 × 106 cells/well in 2 ml of complete culture medium composed of RPMI 1640 medium (Gibco) supplemented with 10% FBS, 1% Penicillin–Streptomycin (Hyclone Laboratories, Inc.) and 1% L‐glutamine (Hyclone Laboratories, Inc.). To mimic a chronic stimulation, cells were cultured for 14 days in a humified CO2 incubator at 37°C in complete culture medium with Dynabeads® Human T‐Activator CD3/CD28 beads (ThermoFisher Scientific) (25 μl/well; bead‐to‐cell ratio of 1:1) and recombinant human interleukin‐2 (rhIL‐2; 25 U/ml; Sigma‐Aldrich). To mimic a transient stimulation, cells were cultured in the presence of CD3/CD28 beads and rhIL‐2 for 3 days and then rested in the presence of rhIL‐2 only for the remaining 11 days of culture. For both stimulation conditions, cells were collected at day 3, 7, 10, and 14. After collection, beads were magnetically removed using a BD IMag Cell Separation Magnet (BD Biosciences) prior to cell preparation for flow cytometry analysis, passaging, cryopreservation. For cell passaging, cells were resuspended in fresh complete medium for chronic or transient stimulation and replated at the same cell density as at day 0.

Flow cytometry

Fresh or thawed T cells were resuspended in 2 ml of BD Pharmingen™ Stain Buffer (FBS; BD Biosciences) and then centrifuged at 300 g for 5 min. For cell surface marker staining, 0.5 − 1 × 106 cells were incubated for 30 min at 37°C with antibody cocktails composed of 100 μl of Stain Buffer (FBS) (BD Biosciences), 10 μl of BD Horizon™ Brilliant Stain Buffer Plus (BSB; BD Biosciences) and each antibody at its recommended concentration, unless otherwise stated. Cells were incubated at room temperature (RT) for 30 min and then washed twice with Stain Buffer (FBS). Cells were resuspended in 0.5 ml of Stain Buffer (FBS) and incubated with the viability dye 7‐AAD (BD Biosciences) at RT for 10 min prior to acquisition on a 3‐laser, 12‐color BD FACSLyric™ Research System or a 5‐laser, 18‐color BD LSRFortessa™ X‐20 Research Use Only system. Instrument configuration details are reported in Tables S5 and S6. For intracellular cytokine detection, the collected cells were first stimulated for 4 h with phorbol 12‐myristate 13‐acetate (PMA; 50 ng/ml; Sigma‐Aldrich) and Ionomycin (500 ng/ml; Sigma‐Aldrich) in the presence of the transporter inhibitors BD GolgiPlug™ and GolgiStop™, as per manufacurer's instructions (BD Biosciences). Cells were then washed with Stain Buffer (FBS) prior to surface marker staining, as per the protocol described above. Cells were then washed in Phosphate Buffered Saline (PBS) without FBS and stained with Fixable Viability Stain 620 (FVS620; BD Biosciences) as per manufacturer's instructions. After two washes in Stain Buffer (FBS), cells were fixed and permeabilized using BD Cytofix/Cytoperm™ Fixation/Permeabilization Solution, as per manufacturer's instructions (BD Biosciences). Cells were then incubated at RT for 30 min with the antibody cocktail composed of 100 μl of permeabilization buffer, 10 μl of BSB and each antibody at recommended concentration, unless otherwise stated. Cells were then washed twice with permeabilization buffer, resuspended in 0.5 ml of Stain Buffer (FBS) and acquired on a 3‐laser, 12‐color BD FACSLyric Research System. After acquisition, all data were exported as FCS3.1 files and analyzed using FlowJo™ software (version 10.6, BD Biosciences). The following mouse anti‐human antibodies, all provided by BD Biosciences, were used in this study: CD4 APC‐H7 (clone RPA‐T4; 5 ul per test), CD4 APC‐R700 (clone SK3; 5 ul per test), CD4 BUV805 (clone SK3; 5 ul per test), CD8 APC‐H7 (clone SK1; 5 ul per test), CD8 AlexaFluor® 700 (clone RPA‐T8; 5 ul per test), CD8 BUV395 (RPA‐T8; 5 ul per test), CD223 BV480 (LAG‐3; clone T47‐530; 0.5 ul per test) CD223 AlexaFluor® 647 (LAG‐3; clone T47‐530; 5 ul per test), CD45RA APC‐H7 (clone HI100; 5 ul per test), CD62L FITC (clone DREG‐56; 20 ul per test), CD95 BV786 (clone DX2; 0.5 ul per test), CD366 BV711 (TIM‐3; clone 7D3; 5 ul per test), CD357 BV421 (GITR; clone V27‐580; 5 ul per test), CD152 PE (CTLA‐4; clone BNI3; 20 ul per test), CD39 BUV737 (clone Tu66; 0.5 ul per test), CD103 APC (clone Ber‐Act8; 5 ul per test), interferon gamma FITC (IFN‐gamma; clone B27; 20uL per test), interleukin‐2 PE (IL‐2; clone MQ1‐17H12; 20 ul per test), and Tumor Necrosis Factor APC (TNF; Mab11; 0.25 ul per test). Biological controls (unstimulated T‐cells) and fluorescence minus one (FMO) controls were used to set gates.

Single cell labeling with sample tags and AbSeq

Cell surface staining was performed as described in the protocol “Single Cell Labelling with the BD Single‐Cell Multiplexing Kit and BD AbSeq Ab‐oligos” (BD Biosciences). Briefly, cryopreserved T cells from three donors, were thawed as per the protocol described in cell processing section. To enable all samples (N = 5 time points) for each donor to be loaded on a single BD Rhapsody™ cartridge, the BD™ Human Single‐Cell Multiplexing kit (BD Biosciences) was used to label the cells from each donor with unique sample tag. Cells were sequentially labeled with sample tags followed by BD AbSeq antibody‐oligos (Ab‐oligos, listed in Table 1; all were used at a volume of 2.5 ul/test, as per the manufacturer's recommendation). First, 1 million cells from each donor/stimulation condition were transferred to a vial containing a unique sample tag barcode (per donor) and incubated at room temperature for 20 min. Following incubation, cells were washed 3 times with Stain Buffer (FBS). Cells were counted and pooled together at an equal ratio of all conditions for each donor. A panel of Ab‐oligos, described in Table 1, was prepared and added to the tube of 1 million pooled cells from each donor. Cells were incubated on ice for 30 min. Following incubation, cells were washed 3 times with Stain Buffer (FBS) and resuspended in Sample Buffer (FBS).

Single cell capture and cDNA synthesis

Cell capture was performed as described in the protocol “Single Cell Capture and cDNA Synthesis with the BD Rhapsody Single‐Cell Analysis System” (BD Biosciences), using both the BD Rhapsody Scanner and Express instrument. The BD Rhapsody scanner was used to perform cell count and viability using Calcein AM (Thermo Fisher Scientific) and Draq7 (BD. Biosciences). 20,000 pooled cells from each donor were loaded into three separate BD Rhapsody cartridges followed by cell capture beads. 13,000 to 16,000 cells were captured per patient, with an average doublet rate of 4.15%. Cells were lysed and the capture beads were then retrieved and washed. Reverse transcription, followed by Exonuclease I treatment was performed on the retrieved cell capture beads, following manufacturer's instructions.

Library preparation

After undergoing cell capture and reverse transcription Rhapsody beads were taken into library preparation as described in “mRNA Targeted, Sample Tag, and BD AbSeq Library Preparation with the BD Rhapsody Targeted mRNA and AbSeq Amplification Kit”. Briefly, all beads were amplified in PCR1 using the BD Human Immune Response Panel (399 amplicons) + SMK, for 11 PCR cycles. Post‐PCR1 reaction cleanup utilized a double‐sided 0.7×/1.2× Ampure method to separate the larger mRNA PCR products from the smaller AbSeq/sample tag PCR products. Each of the Sample tag and mRNA products were taken into separate PCR2 reactions that utilize universal (SMK) or nested (mRNA) primers. 1.2× and 0.8× single‐sided Ampure cleanups, respectively, were performed on the PCR2 products. AbSeq products were taken directly into indexing PCR after PCR1. mRNA PCR2 products (diluted to 1.4–2.7 ng/μl) and AbSeq PCR1/sample tag PCR2 products (diluted to 1.1 ng/μl) were taken into a 6‐cycle indexing PCR reaction. SMK, mRNA, and AbSeq libraries from the same donor were indexed with the same reverse primer, with distinct indexes used between donors. Indexing PCR products for mRNA and AbSeq SMK utilized 0.7× and 0.8× single‐sided Ampure cleanups, respectively. Indexing PCR reactions for mRNA and AbSeq/SMK All libraries were quantified using Agilent High Sensitivity DNA Analysis kits.

Sequencing

All libraries were diluted to 2 nM before pooling for sequencing. Preliminary sequencing for quality assessment was performed on an Illumina NextSeq 500 using a High Output 150 cycle kit with 75 × 75 bp PE reads. Libraries were pooled at a ratio of 1:5:12.5 (sample tag:mRNA:AbSeq) targeting 400 reads/cell from sample tag libraries, 2000 reads/cell from mRNA libraries, and 5000 reads/cell from AbSeq libraries. Full sequencing was done on an Illumina NovaSeq 6000 using an S1 kit with 75 × 75 bp PE reads. For the NovaSeq run the libraries were pooled at a ratio of 1:13:60 targeting an additional 150 reads/cell for sample tag libraries, 2000 reads/cell for mRNA libraries, and 9000 reads/cell from AbSeq libraries. Sequencing metrics are reported in Table S4.

Bioinformatics analysis

FASTQ files were downloaded from Illumina BaseSpace and uploaded onto the Seven Bridges website. Each sample was run separately through the BD Rhapsody Analysis Pipeline using fastqs from the AbSeq panel and Human Immune Response Panel and using the “Single‐Cell Multiplex Kit–Human” multiplexing setting. Output files in csv format were imported into SeqGeq v1.5 software (BD Biosciences) for AbSeq and scRNA‐Seq data analysis. All analyses of molecular cytometry data were performed on cells gated first for quality control measures. These include: (1) exclusion of empty wells or doublets based on a plot of genes expressed versus library size; (2) exclusion of genes never expressed or expressed by all cells by plotting number of cells expressing each gene versus total reads for each gene; and (3) exclusion of genes that do not vary across cells by plotting dispersion versus the number of cells expressing a gene. These measures are described in SeqGeq documentation. Finally, analyses were pregated on CD3+ cells. Data have been deposited at the following link: https://tinyurl.com/cf5bmduf; data are stored in the NCBI BioProject system, under Submission ID SUB9546922 and Bioproject ID PRJNA726263.

Dimensionality reduction

In order to overcome visualization artifacts associated with sparse data, dimensionality reduction was performed using principal component analysis (PCA) guided t‐distributed stochastic neighbor embedding (tSNE). The Opt‐SNE optimized tSNE calculation was used to automatically detect and implement appropriate settings for this machine learning step in analysis [37]. Principal component analysis was performed on highly dispersed gene parameters in combination with extra‐cellular antibody parameters detected via BD's™ AbSeq pipeline. Data were projected on a relative scale (based on bins or channels), in order to account for the vast differences in magnitude for protein versus transcript expression.

Differential expression analysis

Differential expression analysis was performed by pairwise comparisons in volcano plots; illustrating log2 fold change versus adjusts p‐values, also known as “q‐values”, for differentially expressed genes (DEG). Mann–Whitney U‐tests were utilized to estimate the reproducibility of observations in non‐parametric distributions [38]. False Discovery Rate (FDR) adjusted p‐values were appropriate for clusters greater than 200 events in size. Inclusion criteria for DEG: fold‐change values +/−2.0 (up and down regulated, respectively) and q‐values <0.05.

Heat‐maps

Single‐cell heatmap figures were generated by first down sampling populations to a representative number of events. Event columns were then colored on a 0–100% min‐max pseudocolor scale based on relative parameter expression, annotated in descending order.

CONFLICT OF INTEREST

MC, MN, WEL are employed at BD Biosciences. SS and IT were employed by BD Biosciences when the majority of experimental work and manuscript preparation were performed. PC licensed an unrelated invention to BD Biosciences.

AUTHOR CONTRIBUTIONS

Mirko Corselli: Conceptualization (lead); data curation (lead); formal analysis (lead); investigation (lead); methodology (lead); project administration (supporting); visualization (lead); writing – original draft (supporting); writing – review and editing (supporting). Suraj Saksena: Conceptualization (equal); funding acquisition (lead); project administration (lead); writing – review and editing (supporting). Margaret Nakamoto: Methodology (supporting). Woodrow Lomas: Methodology (supporting). Ian Taylor: Data curation (supporting); methodology (supporting); software (supporting); visualization (lead).

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1002/cyto.a.24496.

MIFlowCyt Compliant Items.

Click here for additional data file.

Figure S1. Comparison of protein expression analysis via flow cytometry and AbSeq. (A) Percentage of CD8+ cells expressing CD95, CD103, LAG‐3, CD45RA, CD62L and CTLA‐4 at day 0 (D0), day 3 (D3), day 7 chronic (D7 C) and day 14 chronic (D14 C) of chronic stimulation using either flow cytometry (red trace) or AbSeq (blue trace). (B) Expression of CD95, CD103, LAG‐3, CD45RA, CD62L and CTLA‐4 on CD8+ T cells at day 0 (D0), day 3 (D3), day 7 (D7 C) and day 14 (D14 C) of chronic stimulation. Mean fluorescence intensity (red trace, left y‐axis) and mean molecular count (blue trace, right y‐axis) were used to measure relative antigen expression levels using flow cytometry and AbSeq, respectively. This side‐by‐side analysis was performed using aliquots of cryopreserved cells derived from the same donor (Donor 1), and collected at the indicated time points. Flow cytometry data were down‐sampled in order to analyze the same number of cells as AbSeq.

Click here for additional data file.

Figure S2. Classification of cells based on expression of mRNA only, protein only, or both. Colors represent stimulation conditions.

Click here for additional data file.

Figure S3. Correlation between gene and protein expression at the single cell level. (A) Kinetic expression of genes and corresponding proteins measured as mean molecular count (MMC; red, left y‐axis = protein expression; blue, right y‐axis = gene expression) at day 0 (D0), day 3 (D3), day 7 chronic (D7 C) and day 14 (D14 C) of chronic stimulation. The analysis was performed on three donors. Data are represented as mean ± SD.

Click here for additional data file.

File S1. Supplementary Information.

Click here for additional data file.

Table S1 List of proteins (Ab) and genes upregulated in CD8+CD45RA+CD28+CD27high naïve cells, as compared to CD8+CD45RA+CD28+CD27low cells (fold‐change ≥2, q‐value ≤0.05).

Table S2. List of proteins (Ab) and genes upregulated in CD8+CD45RA+CD28+CD27low cells, as compared to CD8+CD45RA+CD28+CD27high naïve cells (fold‐change ≥2, q‐value ≤0.05)

Table S3. List of shared proteins (Ab) and genes upregulated in CD8+CD45RA+CD28+CD27low and CD8+CD45RA+CD28− TEMRA cells, as compared to naïve cells (fold change ≥2, q ≤ 0.05).

Table S4. Sequencing metrics.

Table S5. BD LSRFortessa X‐20 configuration.

Table S6. BD FACSLyric configuration

Click here for additional data file.