Multiplexed, Sequential Secretion Analysis of the Same Single Cells Reveals Distinct Effector Response Dynamics Dependent on the Initial Basal State

Abstract The effector response of immune cells dictated by an array of secreted proteins is a highly dynamic process, requiring sequential measurement of all relevant proteins from single cells. Herein, a microchip‐based, 10‐plexed, sequential secretion assay on the same single cells and at the scale of ≈5000 single cells measured simultaneously over 4 time points are shown. It is applied to investigating the time course of single human macrophage response to toll‐like receptor 4 (TLR4) ligand lipopolysaccharide (LPS) and reveals four distinct activation modes for different proteins in single cells. Protein secretion dynamics classifies the cells into two major activation states dependent on the basal state of each cell. Single‐cell RNA sequencing performed on the same samples at the matched time points further demonstrates the existence of two major activation states at the transcriptional level, which are enriched for translation versus inflammatory programs, respectively. These results show a cell‐intrinsic heterogeneous response in a phenotypically homogeneous cell population. This work demonstrates the longitudinal tracking of protein secretion signature in thousands of single cells at multiple time points, providing dynamic information to better understand how individual immune cells react to pathogenic challenges over time and how they together constitute a population response.


Introduction
The advent of high throughput single-cell transcriptomic analysis has enabled global unbiased analysis of gene expression in thousands of individual cells and has revolutionized how we study the biological mechanisms at the whole organism level and the complex physiological systems such as the immune system. [1][2][3] However, accessing protein information disease mechanism and drug/vaccine development. [10] Currently, most commonly used tools for multiplexed single-cell secretion analysis includes ELISpot (Fluorospot), [11] intracellular cytokine staining (ICS) flow cytometry with either fluorescence or mass spectrometry (CyTOF) detection. [12,13] Some other multi plexed proteomics profiling methods with lower throughput were also reported, [4,14,15] for example, Herr and co-workers developed single-cell western blotting with ≈11 protein targets detected per cell after multiple fluorophore bleaching-restaining cycles. [15] However, all these high-plex methods could not retain live cells after the assay, which make them impossible to track the same single cells at different time points but also measure a whole panel of protein secretions. [16] And such type of dynamic assay can provide unique information, which inspires deeper insights resolving how immune cells respond and progress upon stimulations. [17][18][19] With the development of microtechnology, a nanowell-based microengraving assay was reported, in which sequential release of cytokines from polyfunctional human T cells can be captured. [20,21] However, the highest degree of multiplexing is four proteins, which is insufficient to dissect the full functional spectrum of heterogeneous immune cells. [10,22] Although an integrative microfluidic-based device was developed to probe single-cell multiplexed input-output dynamic, the cell number to be probed is low (only dozens) and device manipulation is highly complex. [23] To characterize the dynamic, full spectrum secretion information from large quantities of individual cells to understand immune function diversity, it is desirable to develop a method to measure the same single cells over multiple time points (e.g., ≈4 or more) for a panel ≈10 or more protein secretions and such data can be obtained from a large number (e.g., ≈5000) of single cells while minimize the complexity of device handling.
Here, we described a single-cell microchip thatallows for high throughput (≈5000 single cells), multiplexed (≈10 proteins), sequential (4-5 time points) secretion analysis of the same single cells. It was used to profile homogeneous human macrophages and revealed inherently heterogeneous responses over time upon activation with TLR4 ligand lipopolysaccharide. Importantly, we found the stimulated response exhibited two intrinsic states that appear to be associated with the basal function. Single-cell transcriptomic profiles collected at different time points from the samples stimulated in the same way also confirmed the presence of two cellular states with distinct gene expression profiles. Comparing single-cell protein secretion dynamics and transcriptome sequencing data allows for tracking the states of the same cells that confirmed the two states are dependent on initial cell states and differentially regulated by translational and proinflammatory programs.

Single-Cell Secretomic Analysis Microchip with Higher Throughput
The configuration of microchip platform used for dynamic multiplexed single-cell assay was modified from previously reported devices, [24][25][26][27] which was comprised of two components: a high-density antibody barcode patterned glass substrate for surface immunoassay and a nanoliter microtrough array for single cell capture. Notably, we redesigned flow patterning microchip to combine spatial multiplexing (multi antibodies coflow patterning) and spectral multiplexing (multicolor detection) ( Figure 1A) in a much shorter microtroughs (≈0.48 mm) to achieve significantly higher number of single cells to be assayed simultaneously (more than 5000 single cells data can be obtained in one microchip which is around fivefold increase compared to previous work). The flow patterning microchip for antibody immobilization consists of 120 repetitive barcodes, each of which contains 5 stripes in duplicate. The antibody stripes are 30 µm in width and the pitch size of a full barcode is 250 µm. Due to increased flow resistance with much longer channel length, the antibodies were flow patterned with two separated paths to solve this problem. We validated this new flow patterning strategy with both fluorescent-bovine serum albumin (BSA) and recombinant protein sandwich immunoassay to make sure the reproducibility of antibody barcode is adequate for single-cell experiments (FiguresS1-S3, Supporting Information). The polydimethylsiloxane (PDMS) cell capture chip ( Figure 1B) was designed according to the dimension of flow patterned antibody microarray such that it is long enough to contain at least a full set of barcodes, thereby eliminating the need for precise alignment of the antibody barcode slide and the microtrough array PDMS slab. With this design, more than 5000 single cells data (around 30% of total number of microtroughs (n = 18 000), Figure 1C) can be detected in one microchip with minimal sacrifice of parameters to be plexed. For example, up to 15 different proteins can be profiled at the same time if three color detection strategy were employed.

Multiplexed, Sequential Secretion Analysis from the Same Single Macrophages Reveals Heterogeneous Cytokine Secretion Dynamics
One unique feature of our single-cell assay platform is that cells assayed are alive and still isolated in defined locations (specifically for adherent cells). High reproducibility in protein secretion frequency is also validated ( Figure S4, Supporting Information). All these make it possible to accurately and dynamically track the secreted proteins from the same single cells at different time points. Briefly, after measuring protein secretion from single cells over a period of time, we removed the antibody barcode slide that captured the basal secretion profile, and then replaced with a new antibody barcode slide to measure protein secretion from the same single cells for another period of time, during which stimulation reagents can be added, withdrawn, or combined, permitting flexible design of the experiment to perturb cell signaling but keeping track of the same single cells over time. Repeating this process will lead to the measurement of single-cell protein secretion dynamics (Figure 2A). The PDMS microchip for cell capture is oxygen plasma treated for 1 min just before single cell experiment to make its surface hydrophilic to enhance cell adhesion and minimize nonspecific protein adsorption. [28] When changing a new antibody array slide, the PDMS microtrough chip was rinsed three times (including washing and incubation for 2 min in each step) with fresh medium to wash out residual secreted proteins. This also removes detached cells that may dislodge to neighboring microtroughs. Figure 2B shows that 56% of macrophage cells could be retained after the removal and replacing with a new antibody slide. Figure 2C shows a representative single cell that was retained in the same microchannel throughout the entire secretion dynamics experiment and the corresponding secretion patterns.
We applied this platform to investigate the dynamics of U937 derived macrophages in response to TLR4 ligand lipopolysaccharide, which simulates the innate immune response to Gram-negative bacteria. [29,30] The U937 monocyte was differentiated into macrophages by 50 ng mL −1 PMA for 48 h and confirmed by CD14 surface marker staining ( Figure S5, Supporting Information). [31,32] The antibody pairs used in this study (Table S1, Supporting Information) were validated with corresponding recombinant proteins for crosstalk reactivity to ensure technical validity. We also obtained the titration curves, which demonstrated the feasibility of comparing the amount (or concentration) of secreted proteins semiquantitatively ( Figures S6,S7, Supporting Information). After assaying for four time points, a set of data comprising 1752 single cells, each of which has data for a full time course and simultaneous detection of 10 secreted proteins was successfully obtained. It allowed to compare the dynamic change of different secreted proteins of each cell at basal state and with LPS stimulation at different time points ( Figure 2D). For example, TNF was first secreted upon LPS stimulation and the intensity decreased after a 2 h period. Other effector proteins like IL-6 and IL-10 were secreted afterward, which is in agreement with previous reports. [33,34] The dynamic change of each protein in every single cell during the time course can thus be visualized in Figure 2D by connecting their respective detection results at each time point. All single cells (n = 1752) were classified into four subgroups according to their secretion dynamics: on-off, off-on, all-on, and others. The secretion dynamics designated as "on-off" means the protein secretion is active at early time window(s) but stopped in a later period; "off-on" means the secretion was not detectable initially but became active later on; "all-on" means this specific protein was secreted continuously throughout the time course of observation; "other cases" include no secretion or oscillatory secretion patterns. It is noted that the grouping result differs with respect to specific proteins of interest. We noticed that CXCL8 and CCL2 secretions in most cells is "on" during (at least partially) the period of observation, but their dynamics are quite different: the majority of the cells are not secreting CXCL8 during basal state and start turning on CXCL8 secretion after being stimulated. However, the timing for them to turn on CXCL8 secretion is heterogeneous, with comparative numbers of cells turn on CXCL8 secretion between 0-2, 2-4, and after 4 h after being stimulated. Quite differently, CCL2 secretion in the vast majority of cells (71.1%) are kept in the "on" mode independent of being Adv. Sci. 2019, 6, 1801361  responses in immune cells. According to our data, a majority of the cells are classified as "others," meaning not secreting or secreting in oscillatory manner. This is consistent with the previous reports of NF-κB oscillation during transcriptional regulation. However, we also observe a small portion of cells with stable "on-off" or "off-on" secretion patterns, breaking the rule of oscillation. Interestingly, within the population with stable secretion patterns, there are far more cells showing "off-on" pattern versus "on-off" pattern in IL6 secretion (7.1% vs 0.6%) while the trend is reversed in the case of TNF secretion (4.7% vs 13.4%), suggesting cellular basis of a more transient nature of TNF response. We selected 7 proteins being secreted by noticeable number of cells and with apparent secretion alterations during the course of measurement for this analysis. Figure 2D shows the single-cell secretion patterns of 4 proteins and 3 others are shown in Figure S8 in the Supporting Information. Previous study has revealed the stochastic and variable cell switching dynamics of NF-κB in response to TNF activation at single cell level. [19] Herein we show that the transcriptional and functional output of NF-κB activation such as cytokine secretion could exhibit oscillation patterns but the cytokine function outputs are more diverse, which has never been observed previously.

Responses of Macrophages are Intrinsically Heterogeneous and Dependent on the Basal Functional State
By combining all the single-cell data at four different time points (40 proteins parameters) into a unique dynamic, multiplexed single-cell data metric, which cannot be obtained using other methods, we performed unsupervised hierarchical clustering and resolved two distinct major clusters with 1133 and 619 cells, respectively (Figure 3A), indicating the presence of a dynamic heterogeneity within a phenotypically homogeneous macrophage population in response to LPS. Cluster 1 (red label) cells are more active in secretion than cluster 2 (black label) cells not only after stimulation but also in the basal state. We applied a high-dimensional data analysis tool viSNE to visualize multiplexed single-cell data, [35] from which two similar clusters can be obtained too with 943 and 809 cells, respectively ( Figure 3B). Further analysis identified a ≈70% overlap in both clusters using those two clustering methods, indicating the robustness of results regardless of clustering algorithms. The expression of specific proteins in single cells within each cluster (e.g., CCL4, TNF in Figure 3C, and Figure S9 in the Supporting Information) at different time points can be shown in viSNE plots, from which we confirmed the cells more active at basal state are more prone to secreting more proteins upon pathogen stimulation. This phenomenon was not observed previously, highlighting the value of high throughput dynamic detection platform with single cell resolution in resolving how biological systems respond and evolve.
Quantitative analysis of single-cell polyfunctionality (the ability of a cell to cosecrete multiple cytokines and chemokines simultaneously [36] ) was performed and compared between the basal state and the stimulated states, which showed an increase of highly polyfunctional cells upon stimulation (Figure 4A). It is previously reported that polyfunctionality index (PI) serve as an effective parameter to numerically evaluate the degree of polyfunctionality from multidimensional single cell data, permitting more sophisticated statistical analysis. [37] We applied a previously described approach (detailed in Supporting Information) to calculate PI of single cells. [37] Consistent with previous result, we see an increase of PI immediately after stimulation comparing with cells in resting state ( Figure 4B). Interestingly, PI remains constant in the early stage of activation but increased after 4 h being stimulated. This suggests that macro phages increase polyfunctionality immediately (<2 h) after being stimulated but it takes several hours (>4 h) for them to further increase polyfunctionality an enter a more advanced activation state. We then compared polyfunctionality of the same individual cells between different time points and looked for correlative patterns (Figures 4C-E). Interestingly, the resulting averaged polyfunctionality of stimulated single cells showed a linear correlation with the basal state polyfunctionality (R 2 = 0.99, 0.96, 0.98 at 0-2, 2-4, 4-6 h respectively after LPS stimulation). This result means that the higher the numbers (types) of proteins that a cell is secreting in its basal state, the higher numbers (types) of proteins it will likely secrete after being stimulated, and the basal state secretion activity can be taken as an indication of its later activity, suggesting the activation potential of macrophages likely predetermined even before the arrival of external stimulation and is encoded in the basal state intrinsically (Figures 4D,E).

Single-Cell RNA-Seq Confirms Two Major Clusters with Distinct Gene Expression Profiles
We also performed single-cell RNA sequencing analysis on the same samples (both basal and stimulated for different times) in order to compare to single-cell protein secretion data and further investigate the mechanisms underlying the observed heterogeneous states. A massively parallel 3′ mRNA capture and barcoding in droplets was used and the library construction was similar as previously described. [38] The raw sequencing data were processed and analyzed for quality assessment ( Figure  S10, Supporting Information) and the generation of singlecell gene expression matrix. Single-cell transcriptome data of all samples including both basal and activated macrophages of each time point was combined and analyzed with the R package Seurat. Graph-based clustering (resolution = 0.15) identified 3 clusters when all the samples combined, namely, Clusters 0, 1, and 2 (Figures 5A,B). Cluster 2 largely overlaps with the basal or resting macrophages before activation, while Clusters 0 and 1 largely correspond to activated macrophages, indicating the existence of two major distinct activation states in agreement with single-cell protein secretion data. Those two activation states consistently exist in samples of all time points post activation, marking a highly consistent and possibly stable dichotomy in activation states of human macrophages in response to LPS and the top-ranked signature genes defining each cluster were identified ( Figure 5C and Figure S11, Supporting Information). Although mRNA data of activated macrophages supports protein data well, partially confirming our argument that cellular states are possibly predetermined and stable, we notice that the separation of resting macrophages clearly into two clusters is not seen in the transcriptomic profile. In Figures 5A, the observation of three clusters supports the two-state activation model derived from single-cell protein sequential secretion. As expected for whole transcriptome analysis, tSNE separates basal from stimulated cells. The latter however exhibits two clusters, suggesting two major activation states as revealed by single-cell protein data. Figure 5B further confirmed the gradual change of cell states over time and two major states are always present at any given time points. Interestingly, the basal state cells (Cluster 2) already showed "bifurcation" toward two different activation states as defined by Clusters 0 and 1. All these are in agreement with single-cell protein secretion data and support the conclusion that the activation shows two distinct states.
We then asked what gene pathways and biological processes underlie the previously identified two activation states. Gene Ontology (GO) analysis was performed with the R package Gage using the gene lists and p Values generated from the Seurat program as described above, and the marker genes of each of the 3 clusters were used. Surprisingly, the GO terms highly enriched in Cluster 2 involves mRNA processing, which indicates that mRNA processing related genes are highly expressed in the macrophages to prepare them for activation. Furthermore, we found that the GO terms enriched in the two clusters of activated macrophages (Clusters 0 and 1) involve protein translation (ribosomal genes) and inflammatory response (NF-κB pathway, etc.), respectively ( Figure 5D). The anticorrelation relationship between the activity of protein translation and the ability to mount inflammation appears to be mutually exclusive in a single macrophage cell such that each cell has to opt for only one specific state. However, it is unclear if this is a stochastic process or predetermined by the initial   Figure S7 in the Supporting Information). state according to single-cell transcriptome data. The protein secretion dynamics measured on the same single cells suggests the activation states are dependent on the basal state of each macrophage.
We further conducted integrative analysis of single-cell transcriptome and protein secretion data side-by-side, and observed gene-dependent correlations between single-cell RNA and protein profiles. Figure 5E shows 6 immune function genes measured at both transcriptional and protein levels and other genes of interest are shown in Figure S12 in the Supporting Information. The integrated single-cell transcript/protein data can be generally classified into three groups: first, highly correlated genes, including IL6, CCL2, and TNF; second, genes with consistent dynamics but distinct expression levels, including CXCL8 and CCL4; third, genes with consistent expression levels but distinct dynamics. IL10 is an example in this category ( Figure 5E). IL10 mRNA peaks at 2 h, however, IL10 protein secretion reaches maximum between 4 and6 h. This can be explained by the time lag between transcription and translation.
Although the mechanism requires further investigation, our results shows consistency and discrepancy in a gene-specific manner between single-cell mRNA and protein data, which is in agreement with previous reports, [4] highlighting the value to integrate information from multiple omic levels to yield a comprehensive biological picture.

Discussion
We reported a highly multiplexed (≥10) single-cell sequential secretomic profiling platform to evaluate the dynamic immune response from the same single cells at different time points. With this platform, we found human macrophages in response to TLR4 stimulation exhibit two intrinsically distinct states, which were correlated with the basal functional state in each single cell. We also applied single-cell RNA-Seq analysis to confirm the existence of two major states at the transcriptional level and further revealed that the two activation states   are differentially dictated by the expression of protein translation genes and inflammation-related genes, respectively. It has not been possible to obtain this type of comprehensive protein secretion dynamics information of individual cells using any other methods. An even higher degree of multiplexing can be obtained if additional fluorescence channels were used. While this work focused on human macrophages, this approach can be readily applied to other adherent cell types such as epithelial cells, dendritic cells, fibroblasts, etc., making it a versatile platform well suited for much broader applications. [39][40][41] This platform will provide an accessible tool for more in-depth and comprehensive monitoring of cell function and the cognate proteins. It may have the potential to further evolve and be adopted in clinical and pharmaceutical studies to evaluate cellular function heterogeneity or drug responses, [39,[42][43][44] for example, in the immune system or tumor microenvironment.

Experimental Section
Fabrication of Antibody Barcode Array Chips and Microtrough Array Chips: The mold for both the flow patterning PDMS replica and the sub-nanoliter microarray were silicon master etched with deep reactive-ion etching (DRIE). Detailed protocol has been described elsewhere. [25] The PDMS used was RTV615 from Momentive. A and B were in 10:1 ratio.
Cell Culture and Stimulation: Human U937 cell line was purchased from American Type Culture Collection (ATCC) and cultured in RPMI medium 1640 (Gibco; Invitrogen) supplemented with 10% fetal bovine serum (ATCC). The U937 cells were differentiated with 50 ng mL −1 phorbol 12-myristate 13-acetate (PMA) (Fisher) for 48 h, followed by culture in fresh standard medium for 24 h. The cells were harvested with trypsin for single-cell experiments. Cells were challenged with 100 ng mL −1 LPS (Sigma).
Single-Cell Longitudinal Multiplexed Secretion Assay: PDMS microtrough array was first blocked with 3% BSA solution (Sigma) for 2 h and then rinsed with fresh cell medium. Cells were suspended in fresh medium at concentration of 0.2 million cells per mL, followed by the addition of LPS as described previously. The PDMS microtrough array was placed facing upward, and cell culture media was removed until only a thin layer remained on surface. Cell suspension was pipetted (200 µL) onto the microtrough array. After 5 min, the glass slide with antibody barcode was put on top of the PDMS microtrough array with the antibodypatterned side facing the cell-capturing chambers. Then the two parts were clamped tightly with screws using a custom polycarbonate clamping system. Number and locations of cells were confirmed by optical imaging using Nikon eclipse Ti microscope with an automatic microscope stage. Bright field images were obtained. The assembled microchip was placed in a standard 5% CO 2 incubator at 37 °C during the period of cell secretion. After every 2 h, the microchip assembly was dissembled in fresh media and the antibody barcode slide was removed and rinsed with excessive fresh media followed by 2 min of incubation. After repeating this step for 3 times, the cell capture microarray chip was rinsed with fresh media with stimulants added and a new glass slide with antibody microarray was applied on the cell capture microarray to generate a new microchip assemble, followed by imaging and incubation, as described. On the other hand, the glass slide dissembled previously was developed for 1 h at room temperature by introducing a mixture of biotinylated detection antibodies. The detection antibody mixture consisted of the detection antibodies (Table S1, Supporting Information Appendix) at 0.25 mg mL −1 each in 1:200 suspension in 3% BSA. Following this step, the slide was rinsed with 3% BSA solution. The 200 µL of 1:100 suspension APC dye-labeled streptavidin (Biolegend, 5 µg mL −1 ) were added onto glass slide to detect the 635 nm detection antibody group, followed by incubation of 30 min. Following the BSA blocking, the antibody barcode array slide was rinsed with 1× PBS, 0.5× PBS, and DI water sequentially, dried with forced N 2 gas, and then scanned with a four-laser microarray scanner (Molecular Devices; Genepix 4200A) for protein signal detection. Microtrough array images with cell counts were subsequently matched to their protein signals for further data analysis.
Titration Experiment Using Recombinant Proteins in Microfluidic Chips: The titration curves were obtained using recombinant proteins and measured on antibody barcode chips, similar as the ones used for singlecell protein secretion assay (see above). The antibody barcode glass slide was thermally bonded to a PDMS microchannel array slab and then blocked with 3% BSA/PBS solution for 1 h. Recombinant proteins with different concentrations were introduced into different microchannels followed by incubation at 37 °C for 1 h. After that, a cocktail of detection antibodies was added to complete the immune-sandwich assay using the procedure provided by antibody vendors.
Fluorescence Imaging and Analysis: Genepix 4200A scanners (Molecular Devices) were used to obtain scanned fluorescent images. Two channels, 488 (blue) and 635 (red), were used to collect fluorescence signals. The image was analyzed with Genepix Pro software (Molecular Devices) by loading and aligning the microtrough array template followed by extraction of fluorescence intensity values per antibody per microtrough. Fluorescence results were extracted with the image analysis tools in Genepix Pro, and then matched to each of the microtrough array for cell counts as previously extracted from the optical images.
Analysis of Single-Cell Protein Secretion Data: Cell counting was automatically performed by a C++/QT QML software (Isospeak; Isoplexis). Protein signal data were extracted from the multicolor fluorescent images using GenePix Pro 6.1 (Molecular Devices) by aligning a microtrough array template with feature blocks per antibody per microtrough to the protein signal features. Data were extracted using the image analysis tool to gain the mean photon counts per protein signal bar per microtrough and match to the cell counts from the microtrough array. The cell counting and protein signal data were then matched based on their spatial locations. Only the 1-cell wells and their protein signals were used for downstream data analysis. 0-cell wells and their protein signals were used as on-chip controls to provide a measure of local antibody-specific background and were averaged across region on chip. Secretion threshold of a specific factor is defined as mean of the zero-cell wells of its corresponding antibody plus 3 times standard deviation. Values higher than the threshold are taken as "secretion" while values below it are taken as "no secretion" and are changed to 0. The thresholded data were log2 transformed using log2(x+1) before data visualization. Graphpad Prism 7 was used to generate scatter plots and line graphs. Hierarchical clustering and heamaps were performed in R. Polyfunctionality index is calculated using the following equation where n = 6 is the cutoff in the number of functions studied, F i is the frequency (%) of cells performing i functions, and q is a positive number used to modulate the differential weight assignment of each F i . Here q is set to 1 assuming equal weight of each F i . Single-Cell RNA Sequencing: Single-cell RNA Sequencing was performed following the standard DropSeq protocol as described before in detail. [38] Microfluidic device was built following the exact original design file. Beads used in the experiment with oligo synthesize on surface was purchased from ChemGenes. Reagents used were listed in Table S2 in the Supportring Information Appendix. All oligonucleotides used were identical with which was described previously. [38] Sequencing was carried out using Hiseq2000 with 4 samples pooled into one sequencing lane.
Analysis of Single-Cell Transcriptomic Data: Original fastq data was trimmed and transformed into digital expression matrix following the DropSeq data analysis pipeline described in detail elsewhere. Then a filter was applied to get rid of cells expressing fewer than 500 genes, which are likely low-quality cells. Then R package "Seurat" was applied in performing downstream statistical analysis and data visualization using default settings. [45] Marker gene lists and p Values output from "Seurat" was taken as input for gene ontology and pathway analysis, using packages "Gage" and "Pathview" with default settings. [46,47]

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.