Single-cell intracellular epitope and transcript detection reveals signal transduction dynamics

Summary To further our understanding of how biochemical information flows through cells upon external stimulation, we require single-cell multi-omics methods that concurrently map changes in (phospho)protein levels across signaling networks and the associated gene expression profiles. Here, we present quantification of RNA and intracellular epitopes by sequencing (QuRIE-seq), a droplet-based platform for single-cell RNA and intra- and extracellular (phospho)protein quantification through sequencing. We applied QuRIE-seq to quantify cell-state changes at both the signaling and the transcriptome level after 2-, 4-, 6-, 60-, and 180-min stimulation of the B cell receptor pathway in Burkitt lymphoma cells. Using the multi-omics factor analysis (MOFA+) framework, we delineated changes in single-cell (phospho)protein and gene expression patterns over multiple timescales and revealed the effect of an inhibitory drug (ibrutinib) on signaling and gene expression landscapes.

In brief Rivello et al. develop QuRIE-seq for highthroughput quantification of single-cell mRNA levels in combination with extraand intracellular protein levels via immunostaining with 80 DNA-conjugated antibodies. This method shows the cellular response over time to external perturbations at the surface marker, (phospho)proteome, or gene expression level.

INTRODUCTION
The processing of information from the outside environment through signaling pathways is one of the most fundamental processes determining cellular phenotype, function, and fate. Single-cell multi-modal omics tools are changing our understanding of biology, but we still lack high-throughput single-cell techniques that capture both changes in phosphorylation levels and the much slower ensuing changes in gene expression patterns (Zhu et al., 2020). Targeted methods like REAP-seq (Peterson et al., 2017), CITE-seq (Stoeckius et al., 2017), or ECCITEseq (Mimitou et al., 2019) measure simultaneously transcripts and proteins in single cells by using DNA-tagged antibodies, but are limited to the detection of surface epitopes. Our recent work showed a strategy for plate-based single-cell transcriptome and intracellular protein measurements for six epitopes (Gerlach et al., 2019). Here, we present quantification of RNA and intracellular epitopes by sequencing (QuRIE-seq), a highthroughput droplet-based platform to simultaneously quantify 80 intra-and extracellular (phospho)proteins and the transcriptome from thousands of individual cells.

QuRIE-seq workflow
We developed QuRIE-seq to study signal transduction and downstream transcriptional changes by using a validated panel of DNA-barcoded antibodies targeting components of the B lymphocyte (B cell) antigen receptor (BCR) signaling pathway, as well as cell-cycle and surface markers (Table S1). To allow intracellular epitope detection, we cross-linked cells by using a mix of reversible fixatives (dithiobis(succinimidyl propionate)/ succinimidyl 3-(2-pyridyldithio)propionate [DSP/SPDP]), permeabilized, stained with the DNA-barcoded antibody panel, and co-compartmentalized in nanoliter droplets with single-cell barcoded primer-loaded gel beads (based on the inDrop protocol) (Klein et al., 2015;Zilionis et al., 2017), before further library preparation and sequencing ( Figures 1A-1E). Reversible fixation MOTIVATION External stimuli lead to a range of dynamic intracellular responses. There is a dearth of methods that simultaneously capture the dynamics of this response at the levels of signaling and gene expression networks. We combined high-throughput single-cell RNA sequencing with intracellular epitope detection via DNA-barcoded antibodies, thus opening a window on the dynamics of signal transduction.
of Burkitt lymphoma (BJAB) cells with DSP/SPDP and permeabilization showed similar or higher (phospho)protein signals and signal-to-noise ratios by flow cytometry and were equivalent to conventional paraformaldehyde fixation methods. For this BJAB cancer cell line, the fixation and permeabilization method resulted in a gene-detection rate comparable to that of unfixed single-cell analysis, similar to previously published results in primary keratinocytes (Gerlach et al., 2019). Technical quality checks of the transcriptomic and proteomic data generated with QuRIE-seq for fixed BJAB cells are presented in Figures S1A-S1D. Analysis of primary B cells isolated from peripheral blood mononuclear cells resulted in a good protein library (high protein counts and apparent detection of >70 proteins), but low-quality RNA library ($100 genes/cell and high percentage of mitochondrial counts), prompting us not to analyze this dataset further ( Figures S1E-1H).
BJAB cells are characterized by functional BCR signaling that can be further stimulated by using a polyclonal anti-immunoglobulin antibody (aIg) (Schamel and Reth, 2000), as confirmed by aIg-induced phosphorylation of components of the BCR signal transduction pathway and aIg concentration-dependent secretion of CCL3 chemokine and interleukin-10 (IL-10) and IL-6 cytokines (data not shown). We set out to explore molecular changes induced by aIg treatment of BJAB cells at combined Figure 1. Multi-modal signal transduction in BJAB cells at the single-cell resolution by using quantification of RNA and intracellular epitopes (QuRIE-seq) combined with multi-omics factor analysis (MOFA+) (A-G) Schematic overview of the QuRIE-seq workflow for combined single-cell RNA sequencing and (phospho)protein measurement. (A) The dynamic response of BJAB cells to stimulation for different durations with polyclonal anti-immunoglobulin antibody (aIg). The scheme of part of the targeted proteins in the B cell antigen receptor (BCR) signaling pathway is shown. The timescale of (phospho)proteome changes falls into a minute regime, whereas most significant transcriptome changes are observed after hours. (B) Cells are reversibly cross-linked with DSP/SPDP (fixatives) and permeabilized. (C) DNA-tagged antibodies targeting membrane and intracellular (phospho)proteins are used for immunostaining. (D) Single cells are co-encapsulated in a high-throughput manner in nanoliter droplets with barcoding hydrogel beads. (E) Preparation of two libraries (proteome and transcriptome) is followed by sequencing. (F) Data analysis overview: MOFA+ takes sequencing data matrices from both modalities (transcriptome and proteome) and decomposes these into matrices containing factors that can be analyzed further downstream. (G) UMAP embedding based on the MOFA+ factors (factor 1 to 7) for 4,754 BJAB cells stimulated with aIg for different durations: 0, 2, 4, 6, 60, and 180 min.
2 Cell Reports Methods 1, 100070, September 27, 2021 Report ll OPEN ACCESS (phospho)protein and transcriptome levels across a minute-tohour timescale by applying QuRIE-seq on cells stimulated with aIg for 0, 2, 4, 6, 60, and 180 min. Data pre-processing included stringent quality control analysis and normalization, resulting in a multi-modal dataset of 4,754 cells with matched gene and (phospho)protein expression levels. Finally, the data were scaled and corrected for sequencing depth, mitochondrial content (RNA dataset), and histone H3 levels (protein dataset).
Multi-omics factor analysis Next, we used unbiased Multi-Omics Factor Analysis (MOFA+) (Argelaguet et al., 2020) to combine the quantification of (phospho)proteins and transcripts into a single model that identifies factors explaining variation from both modalities over the full time series ( Figure 1F and S2A-S2C). First, we used the computed MOFA+ factors as input for uniform manifold approximation and projection (UMAP), capturing the cellular responses to induced activation of the BCR pathway. As shown in Figure 1G, two distinct phases of cellular responses in BJAB cells are observed after short (0-2 min) and long (60-180 min) stimulation. Principal-component analysis (PCA) of the separate datasets, followed by integration of the modalities by using weighted nearest-neighbor (WNN) analysis from the Seurat v.4 package (Hao et al., 2021), confirmed these observations, showing changes at the minute timescale for the (phospho)protein dataset and the hour timescale for the transcriptome dataset in the first principal component ( Figures S2D-S2H), as expected. Together, these initial observations indicate that the computed factors of QuRIE-seq measurements capture the cellular responses to a signal across modalities and timescales.

Revealing signal transduction dynamics
To explore which molecular changes underlie the BJAB cells' response to the aIg stimulus at the minute and hour timescales, we analyzed the computed MOFA+ factors in more detail and found that factors 1 and 3 most strongly correlated with the time of the treatment ( Figure S2C). Factor 1 captures minute timescale changes ( Figure 2A), with predominant effects on the phosphoprotein levels ( Figure 2B), as confirmed by flow cytometry. Only a minor part of factor 1 is associated with the variance in the mRNA dataset, as transcriptional changes after external stimulation typically take place at the hour, not minute, timescale ( Figure S2A).
Examination of the (phospho)proteins and genes with the highest loading values for factor 1 highlighted components of the BCR signaling pathway (e.g., p-CD79a, p-PLC-g2, p-SRC, and p-SYK, IGHM, IGKC), indicating that factor 1 most likely represents the core elements of the BCR signal transduction pathway, which is activated within 2 min of BCR stimulation with aIg (Figures 2C and 2D). Surprisingly, using our unbiased approach, we identified p-JAK1 as strongly contributing to factor 1, increasing within the first 2 min after aIg addition, suggesting regulation of JAK1 phosphorylation by BCR activation (Figures 2B and 2D). Although the tyrosine kinase JAK1 has been suggested as a substrate for SYK (Patterson et al., 2015;Satpathy et al., 2015), phosphorylation of JAK1 upon BCR activation has thus far not been described. The unexpected activation of JAK1 by aIg stimulation was further confirmed with flow cytometry ( Figure S3A).
Notably, canonical activation of JAK1 by indirect activation of signaling through cytokines is unlikely to explain these effects, as the observed timescale of 2 min is too short for expression and/or secretion of endocrine cytokines. In addition, we found that the JAK1/2 inhibitor ruxolitinib resulted in decreased phosphorylation of the BCR pathway components SYK, Bruton's tyrosine kinase (BTK), and PLC-g2, as well as the JAK1 substrate STAT6 ( Figure S3B). This suggests a previously unidentified role for JAK1 in the BCR response to aIg stimulation revealed by QuRIE-seq. We performed additional analyses on p-JAK high versus p-JAK low cells (top and bottom 5% at each time point). At several time points, we found a common set of phosphoproteins that are differentially expressed in the high versus low cell subpopulations, including p-CD79a, p-SYK, and p-SRC, all related to BCR signaling activation. Even though no corresponding consistent differences in gene expression were found, these findings illustrate that we can distinguish cells with higher versus lower signaling responses to the stimulus within the population of cells ( Figure S4).
Whereas factor 1 captures changes at the minute timescale, factor 3 distinguishes samples at the hour timescale (60 and 180 min, Figure 2E). The (phospho)proteins with the highest factor 3 loading values are p-S6, XBP1, p-p38, p-AKT, p-TOR, and p-STAT5 ( Figure 2F). These (phospho)proteins have been shown to contribute to the increase in mis/unfolded proteins in activated B cells through the unfolded protein response (UPR) (Gaudette et al., 2020), activation of glycolytic energy metabolism, and cell-cycle entry. RNA features with positive factor 3 loading values ( Figures 2F and 2G) are significantly enriched for gene sets associated with similar cellular functions: proteasomemediated protein catabolic process, COPII-mediated vesicle budding, Golgi vesicle organization, and cell-cycle entry. These biological processes relate to increased protein production and energy metabolism, antibody secretion (plasma cell differentiation), and cell-cycle entry ( Figure 2G). In addition to proteasome-mediated protein catabolic processes, the protein folding gene ontology (GO) term is enriched (adjusted p = 0.021) in the RNA dataset, consistent with the enrichment of UPR-related (phospho)proteins associated with factor 3. The expression of genes with the highest factor 3 loading, including NEAT1, NPM1, and BTF3, is upregulated after 60 and 180 min of BCR signaling stimulation with aIg ( Figure 2H). Comparison of differentially expressed genes determined by bulk RNA sequencing and by QuRIE-seq shows similar upregulation profiles at 60 or 180 min (data not shown). In summary, the molecular changes at the hour timescale seem to be dominated by proteins and genes involved in processes related to increased protein synthesis, activated glycolysis, B cell activation, and proliferation.

Disruption of signal transduction by using ibrutinib
Finally, as QuRIE-seq allows us to map signal transduction at both the (phospho)protein and the mRNA level, we explored whether the developed platform can be used to obtain a highresolution view of the mechanism of action of drugs. Therefore, we stimulated BJAB cells with aIg, as described above, in the presence or absence of 1 mM ibrutinib (Ibru), a small-molecule inhibitor of BTK, blocking a part of the BCR signaling ( Figure S5A). We explored the effect of Ibru treatment at both minute and hour  Figure 3C) suggest these processes are, at least partially, BTK dependent. Indeed, phosphorylation of PLC-g2 and BLNK is inhibited by Ibru, as were CD79a, SYK, and JAK1 ( Figure S5B). The decreased phosphorylation of proteins upstream of BTK suggests the existence of a feedback loop that is affected by Ibru treatment. The expression of the top three genes with high factor 3 loading shows diminished upregulation after 180 min of treatment with aIg ( Figure S5C). To identify additional effects of BTK inhibition, we note that factor 5 most strongly correlates (Pearson correlation = 0.47) with Ibru treatment and shows increased values after 180 min of stimulation with aIg for Ibru-inhibited cells compared with cells without Ibru inhibition ( Figures  3B and 3C).
GO-term analysis shows that features highly contributing to this factor in the positive loadings are two negative regulators of G-protein-coupled signaling, RGS2 and RGS13, for the RNA dataset (GO: 0045744, adjusted p = 0.22). RGS2 and RGS13 expression after disturbed BCR activation might relate to the germinal center (GC) phenotype of BJAB cells, where signaling through the G-protein-coupled receptors CXCR4 and CXCR5 orchestrates GC dynamics (Wu et al., 2019). Furthermore, in the protein dataset, p-ERK 1/2 contributes, albeit modestly, to factor 5 ( Figures 3D-3F). Surprisingly, ERK 1/2 phosphorylation is maintained independent of Ibru-mediated BTK inhibition (Figure 3F), which supports the notion that Ibru only partially blocks B cell signal transduction. In line with this, Ibru does indeed inhibit, in a dose-dependent manner, secretion of IL-10 and CCL3, but not IL-6 ( Figure S5D), supporting the notion of BTKdependent and -independent aIg-induced activation of BJAB cells. These precursory findings illustrate the potential of QuRIE-seq to study the complexity of inhibitory drug effects on signal transduction.

DISCUSSION
In conclusion, QuRIE-seq combined with a time series of singlecell intracellular epitope and transcriptome data provide a powerful method to study multi-modal intracellular signal transduction at multiple timescales and characterize the mechanism of action of drugs on signal transduction. Here, we used an antibody panel targeting 80 (phospho)proteins, but we believe there are no biophysical or biochemical limitations to expanding this panel to include 500+ target epitopes. Targeting a large number of (phospho)proteins can reveal previously unidentified signaling networks, such as the involvement of JAK1 in BCR signaling or the BTK-independent activation of ERK 1/2. Although we have exploited droplet-based microfluidics to implement QuRIEseq, we envision future adaptation to high-throughput platebased methods such as sci-Plex (Zheng et al., 2017) and Seq-Well (Gierahn et al., 2017). The MOFA+ framework allowed us to analyze both modalities simultaneously, with clear and interpretable results. With an increasing number of multi-modal single-cell technologies becoming available, complementary analysis methods are rapidly being developed, and other options will be available for analyzing our QuRIE-seq data. For example, we tested WNN (implemented in Seurat v.4;Hao et al., 2021) analysis of our data, and combined PCA and WNN analysis shows results very similar to those of MOFA+ in terms of resolving the time-dependent manner of phosphoprotein activation at short timescales, as well as transcriptional changes at later time points. The combination of customized antibody panels, the ability to detect changes with a high temporal resolution, and the multi-omic readout of the effect of inhibitory drugs leads us to anticipate a broad implementation of QuRIE-seq in fundamental signaling studies and drug development research.
Limitations of the study QuRIE-seq relies on binding of DNA-tagged antibodies to membrane and intracellular epitopes. Cell cross-linking might change protein conformation and thus affect the binding of antibodies, limiting our ability to quantify protein levels of those targets. In this study, we used some 80 antibodies in our panel. Expansion of the panel brings significant additional costs.
Our initial studies on primary B cells indicated that mRNA extraction was insufficient to obtain a good-quality mRNA library. Further improvements in the reversible cross-linking chemistry are required.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following: Antibodies used for immunostaining of cells prior to microfluidic encapsulation are listed in Table S1 See Table S1 N/A Antibodies used for flow cytometry analysis of BJABs are listed in Table S2 See

RESOURCE AVAILABILITY
Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Prof. Dr. Wilhelm T. S. Huck (w.huck@science.ru.nl).

Materials availability
This study did not generate new unique reagents.
Data and code availability d Single-cell RNA-seq data have been deposited at GEO and are publicly available as of the date of publication. Accession numbers are listed in the key resources table. Flow cytometry data reported in this paper will be shared by the lead contact upon request. d All original code has been deposited at Zenodo and is publicly available as of the date of publication. DOIs are listed in the key resources table. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL
The BJAB cell line used in this study were obtained from DSMZ, Germany. The cell line was established in 1973 from the inguinal tumor of an EBV-negative 5-year-old African girl with Burkitt lymphoma at the preterminal stage; although the tumor tissue was experimentally exposed to EBV, neither EBV sequences nor EBV antigens were detectable in the established cell line. Exome and RNA sequence data are available DSMZ (https://www.dsmz.de/collection/catalogue/details/culture/ACC-757).

Antibody labeling
All antibodies (Abs) were purchased purified (see Table S1 for a full list of the Abs used). Abs used in the QuRIE-seq experiment were first validated using flow cytometry.
QuRIE-seq requires the use of DNA-tagged Abs for staining cells prior to encapsulation to characterize the single-cell proteome. In the following is briefly detailed the DNA conjugation strategy to the Abs which made us of standard NHS chemistry. First, Abs were buffer exchanged into 0.2M NaHCO 3 (pH 8.3) (Sigma Aldrich, USA) using Zeba TM Spin Desalting Columns, 40K MWCO (Thermo Fisher Scientific, USA). Then, dibenzocyclooctyne-S-S-N-hydroxysuccinimidyl ester (Sigma Aldrich, USA) was resuspended in DMSO (Sigma Aldrich, USA), added in a 10x molar excess to the Abs and incubated for 2 hours at room temperature. Abs were then buffer exchanged into PBS (Thermo Fisher Scientific, USA) and excess linker molecules were removed using 30K Amicon centrifuge units (Merck, USA). Subsequently, 5' Azide modified oligos (Biolegio, The Netherlands) were added in a 3x molar excess to the chemically modified Abs and incubated for sixteen hours at 4 C in dark. After incubation, Abs were buffer exchanged into PBS (Thermo Fisher Scientific, USA) containing 0.05% sodium azide (Sigma Aldrich, USA) and 0.1mM EDTA (Sigma Aldrich, USA), and excess oligos were removed using 100K Amicon centrifuge filters (Merck, USA). After labelling, DNA-tagged Abs were stored at 4 C. Labelling efficiency was determined by non-reducing SDS page gel analysis where 0.5 mg antibody was diluted in 1x Laemmli sample buffer (Bio-Rad, USA) and loaded on a mini-PROTEAN TGX stain-free gel (Bio-Rad, USA). Gels were imaged on a ChemiDoc Touch imaging system (Bio-Rad, USA).

BJAB cells preparation
BJAB cells were obtained from DSMZ, Germany (ACC 757). Cells were cultured in RPMI medium (Thermo Fisher Scientific, USA) supplemented with 20% FBS (Biowest, France) and 1% Penicillin-Streptomycin (Thermo Fisher Scientific, USA). Before stimulation, cells were harvested, centrifuged at 200 g for 5 minutes and re-suspended in complete medium at a concentration of 3*10 6 cells/ml. For each time point, 500 ml cell suspension was dispensed in 3 ml FACS tubes (BD biosciences, USA). For cells treated with inhibitor 1 mM Ibrutinib (Selleck Chemicals, USA) was added directly to the cell suspension. Subsequently, cells were allowed to rest in the complete medium for one hour at 37 C with or without an inhibitor. Cells were then stimulated with 10 mg/ml F(ab')₂ Fragment Goat Anti-Human IgA + IgG + IgM (H+L) (Jackson ImmunoResearch, USA) for the indicated time points. After incubation, cells were immediately fixated by adding two times concentrated fixative (5 mM dithiobis(succinimidyl propionate) (DSP) (Thermo Fisher Scientific, USA) and 5 mM succinimidyl 3-(2-pyridyldithio)propionate (SPDP) (Thermo Fisher Scientific, USA) in PBS (Thermo Fisher Scientific, USA), and incubated for one hour at room temperature. Fixated cells were quenched and permeabilized with 100 mM Tris-HCl pH 7.5, 150 mM NaCl and 0.1% Triton X100 (Thermo Fisher Scientific, USA) for 10 minutes at room temperature. After permeabilization, cells were washed once and blocked for 45 minutes in Cell Reports Methods 1, 100070, September 27, 2021 e2 Report ll OPEN ACCESS 0.5X protein-free blocking buffer (Thermo Fisher Scientific, USA) with 0.2 mg/ml dextran sulfate (Sigma Aldrich, USA) and 0.5 U/ ml RNasin plus (Promega, USA) in PBS (Thermo Fisher Scientific, USA). Cells were stained in a blocking buffer containing DNAtagged antibodies (Ab list used for QuRIE-seq encapsulation in Table S1, Abs used for flow cytometry in Table S2) for one hour at room temperature. Following staining, cells were washed twice with the blocking buffer and re-suspended in PBS containing 5 mg/ml BSA (Thermo Fisher Scientific, USA) and 0.5 U/ml RNsing plus (Promega, USA). Cells were sorted using the BD Biosciences FACS Melody, based on singlets, into protein LoBind tubes (Eppendorf, Germany). Just before encapsulation, each cell suspension was brought to 10 5 cells/ml and supplemented with 15% (v/v) Optiprep (Sigma Aldrich, USA).

Primary B-cells preparation
Buffy coats (Sanquin, The Netherlands) were firstly divided into two 50 ml tubes, then 750ml RosetteSep B-cell enrichment cocktail (Stemcell, Canada) was added, and blood was incubated for 20 min on a roller bench at room temperature. Subsequently, B-cell enriched PBMCs were isolated by Ficoll density gradient separation (GE Healthcare, USA). Enriched PBMCs were washed twice in PBS (Thermo Fisher Scientific, USA) containing 2mM EDTA (Lonza, Switzerland). Total B-cells were further isolated by the B-cell isolation kit II (MACS, Miltenyi Biotec, Germany) according to manufactures instructions. Purified B-cells were resuspended at a concentration of 4*10 6 cells/ml in RPMI (Thermo Fisher Scientific, USA) supplemented with 1% Penicillin-Streptomycin (Thermo Fisher Scientific, USA) and 50 mM 2-Mercaptoethanol (Thermo Fisher Scientific, USA). 200ml of cell suspension was pipetted into FACS tubes (BD Bioscience, USA) and cells were allowed to rest for one hour at 37 C. Subsequently, cells were stimulated with 10 mg/ml F(ab')₂ Fragment Goat Anti-Human IgA + IgG + IgM (H+L) (Jackson ImmunoResearch, USA) for the indicated time points. After incubation, cells were immediately fixed in an equal volume 2x concentrated fixative, permeabilized and labelled as described for the BJAB cells above.
Flow cytometry analysis DSP-and SPDP-fixed and permeabilized BJAB cells were centrifuged at 800 g and resuspended in PBS (Thermo Fisher Scientific, USA) supplemented with 0.1% BSA (Thermo Fisher Scientific, USA) containing the indicated Abs, and incubated for one hour at room temperature in dark (Table S2 shows different Abs used for flow cytometry analysis). After incubation, cells were washed three times with PBS containing 0.1% BSA and resuspended in PBS with 0.1% BSA. Finally, cells were measured on a FACSverse (BD Bioscience, USA). Flow cytometry results were analysed using FlowJo software (BD Bioscience, USA) (10000 cells we measured per sample).
ELISA BJAB cells were cultured in round bottom 96-wells plates at a concentration of 10 5 cells/well in a volume of 200 ml. Cells were stimulated with 0-20 mg/ml F(ab') ₂ Fragment Goat Anti-Human IgA + IgG + IgM (H+L) (Jackson ImmunoResearch, USA) and incubated for 16 hours at 37 C. The supernatant was collected and analysed for the presence of IL6, IL10 and CCL3 using uncoated ELISA kit (Thermo Fisher Scientific, USA). ELISA was performed according to the manufacturer's instructions.
Bulk mRNA sequencing BJAB cells were harvested, centrifuged and resuspended in complete medium at 0.5*10 6 cells/ml. Thereafter, 4 ml of cells were seeded in T25 flasks (Thermo Fisher Scientific, USA), and cells were incubated for 1 hour at 37 C. After incubation, the cells were stimulated with 10 mg/ml F(ab')₂ Fragment Goat Anti-Human IgA + IgG + IgM (H+L) (Jackson ImmunoResearch, USA) for the indicated time points. After stimulation, cells were immediately collected on ice and centrifuged at 200 g for 5 minutes at 4 C. Supernatant was collected and stored at -20 C for further analysis and the cell pellet was resuspended in the lysis buffer RLT (Qiagen, The Netherlands). Cell lysates were stored at -80 C until further use. Library preparation and sequencing was performed by Macrogen (Macrogen, Korea).

Microfluidic device fabrication for single-cell barcoding
Microfluidic devices were fabricated by using photo-and soft lithography based on the AutoCAD files with channel designs that were kindly provided by Dr. Linas Mazutis (Klein et al., 2015). Silicon wafers were spin-coated with a uniform layer of SU8-2050 photoresist (MicroChem Co., USA), soft-baked, UV exposed through transparency mask (JD Phototools, UK), baked again post-exposure, developed and hard-baked according to manufacturer's protocol (MicroChem Co., USA). After production, the height of the structures was measured by Dektak profilometer (the average height $70 mm). The obtained wafers were used as masters for the PDMS devices. PDMS prepolymer and crosslinking agent were mixed at 10:1 ratio (w/w) and poured on the master bearing the microchannel structure, degassed for 30 minutes and cured at 65 C for at least 2 hours. Subsequently, the PDMS replica was cut out from the master and the holes were punched at inlet and outlet ports using a biopsy puncher of 1 mm inner diameter (pfm medical, USA). Thereafter, glass slides were carefully washed with a mixture of soap and water, and then both the replicas and the glass slides were rinsed with ethanol. The clean replicas were bonded to the glass slides after oxygen plasma treatment (Femto, Diener electronic, Germany). The channels of the devices were rendered hydrophobic by silanization with 2% 1H,1H,2H,2H-perfluoro-octlytriethoxysilane (Sigma-Aldrich, USA) in FC-40 (Sigma Aldrich, USA). The devices were incubated at 95 C overnight and afterward stored at room temperature until further usage for droplet production. Microfluidic single-cell barcoding The chip design described by Klein et al. 6 for single-cell barcoding was used, with four inlets (reverse transcription mix, beads, cell suspension and oil phase with surfactant). The microfluidic single-cell barcoding experimental workflow was followed as described by (Zilionis et al., 2017) with some modifications. Prior to microfluidic encapsulation, the dispersed and continuous phases used for emulsification were prepared. Barcoded hydrogel beads (1CellBio, USA) were spun down, and mixed with 2X concentration bead mix (0.2% vol/vol Triton X-100 and 2x First Strand buffer) and centrifuged at 5000 g for 2 minutes. The supernatant was removed and the concentrated beads were loaded into the tubing with a previously made air spacer to avoid spreading of beads in HFE-7500 fluorinated fluid (3M) loaded syringe. The cell suspension mixed with 15% (v/v) Optiprep (Sigma-Aldrich, USA) was loaded using a pipette-tip based method for seeding cells to droplet microfluidic platforms, described elsewhere (Sinha et al., 2019). The reverse transcription (RT) mix comprising 1.3X RT premix MgCl 2 (1M), DTT (1M), RNAse OUT (40U/mL), Superscript III (200U/mL) and RNAse free water (Ambion, Thermo Fisher Scientific, USA) was loaded into the tubing. The tubing with the RT mix was wrapped around the ice pack to prevent the increase of temperature and deactivation of enzymes. The oil phase consisted of 3x (v/v) diluted 10% (w/w) RAN (RAN Biotechnologies, USA) in HFE-7500 oil. Solutions were introduced into the microfluidic chip using neMESYS (Cetoni GmbH, Germany) syringe pumps. First, the beads were injected at $200 mL/h. As soon as the beads appeared at the inlet port, the flow was decreased down to $20 mL/h, and cell suspension, RT mix, oil spacer were loaded at flow rates $150 mL/h. When all reagent started to flow through the channels, the final flow rates of the beads, cell suspension, RT mix, oil phase were established $15-20 mL/h (the flow was tuned during the experiment), 100 mL/h, 100 mL/h, and 90 mL/h, respectively, to produce 2 nl drops. Closely packed barcoding beads were co-encapsulated with single cells in nanolitre droplets, reaching 60-90% cell barcoding efficiency.
After microfluidic emulsification, the cells were de-crosslinked by the action of DTT present in the RT mix which reduced the disulphide bonds present and the barcoding primers were released from the hydrogel beads by 8-minute UV exposure at 10 mW/ cm2 from a photolithography machine (ABN Inc, USA). Subsequently, the poly(dT) tail of the barcoding primers annealed both to the poly(A) tail of the mRNA and the poly(dA) tail of the Ab tags and the single cells-derived material was barcoded by reverse transcription: incubation of 2 hours at 50 C, followed by 15 minutes at 70 C. The emulsion was then broken by adding 3 mL of 20% (v/v) PFO on the top of the emulsion, and then 40 mL of HFE-7500 fluorinated fluid (3M). The samples were stored until further use at -80 C.

Size separation of mRNA and protein libraries
After microfluidic single-cell barcoding, samples were thawed, centrifuged at 4 C for 5 minutes at 19000 g to pellet cell debris, then the aqueous post-RT material (top layer) was separated from the bottom oil layer and transferred to a nucleic acid purification column (Corning Costar Spin-X column, 0.45 mm, Sigma-Aldrich, USA). Tubes were centrifuged for 1 minute at 16000 g at 4 C. Subsequently, 100 mL of digestion mix was added per each 70 mL of the sample (final concentrations were 0.5x FastDigest Buffer (Thermo Fisher Scientific, USA) and 0.6 U/mL ExoI (Thermo Fisher Scientific, USA)) and incubated for 30 minutes at 37 C. After digestion of unused primers and primer dimers, the mix was size separated using 0.6x volume of AMPure XP magnetic beads (Beckman Coulter, USA). The short fragments (the proteomic library that was retained in the supernatant) were separated from the long fragments (the transcriptomic library kept by beads) and each library was then processed independently.

mRNA library preparation
After the size separation of the mRNA library from the protein library, the mRNA library was processed for second-strand synthesis (SSS) and linear amplification by In Vitro Transcription (IVT). Briefly, 2 mL of 10x SSS buffer and 1 mL of SSS Enzyme (NEBNext mRNA Second Strand Synthesis Module, New England Biolabs, USA) were added to the 17 mL of the purified sample, mixed, and incubated at 16 C for 1.5 hour, and 20 min at 65 C. Subsequently, 60 mL of IVT mix was added to each of the 20 mL SSS product, and incubated at 37 C for 15 hours. The IVT reaction was composed by 1x reaction buffer, 10mM each of ATP, CTP, GTP and UTP, and 10x diluted T7 RNA Polymerase mix (HiScribe TM T7 High Yield RNA Synthesis Kit, New England Biolabs, USA). After linear amplification, the amplified RNA material was purified with 1.3x AMPure XP magnetic beads (Beckman Coulter, USA), and fragmented by adding 1 mL of 10X RNA fragmentation Reagents (Thermo Fisher Scientific, USA) to 10ul of purified IVT, mixing and immediately incubating at 70 C for exactly 3 minutes. Thereafter, the mix was transferred on ice and 34 ml of STOP mix (8.9 ml nuclease-free water, 24 ml AMPure XP magnetic beads, and 1.1 ml STOP solution) were added to it. Purification using AMPure XP magnetic bead was followed by sample elution in 8 ml nuclease-free water. The following step of the library preparation consisted of the reverse transcription (RT) of the amplified RNA material using random hexamers. This was performed by first addition of 1 mL of 10 nM of each dNTP (Thermo Fisher Scientific, USA) and 2 mL of 100 mM of PE2-N6 primer (Biolegio, The Netherlands), and incubation of the mix at 70 C for 3 minutes, followed by cooling down the mixture on ice. Next, 3.5 mL nuclease-free water was added to 4 mL of 5x PrimeScript buffer (Takara, USA), 1 mL RNAse OUT (Thermo Fisher Scientific, USA) and 0.5 mL PrimeScript RT (Takara, USA), and the reaction mix was incubated at 30 C for 10 minutes, followed by 42 C for 1 hour and by 15 minutes at 70 C. The RT product was then purified with 1.2x AMPure XP magnetic beads and eluted in 10 mL of nuclease-free water and amplified by PCR. The number of required cycles was determined by qPCR using 0.5 mL of purified material with the addition of 6.5 mL nuclease-free water, 10 mL of 2x Kapa HiFi HotStart PCR mix (Roche, Switzerland), 1 mL of 20x EvaGreen Dye (Thermo Fisher Scientific, USA), and 2 mL of PE1/PE2 primer mix (5 mM each) (Biolegio, The Netherlands). Final material was obtained by mixing 0.5 mL of nuclease-free water, 12.5 mL of 2x Kapa HiFi HotStart PCR mix (Roche, Switzerland), and 2 mL of PE1/PE2 primer mix (5 mM each) (Biolegio, The Netherlands) and the 9.5 mL of the sample. The PCR thermal cycling was as follows: 2 minutes at 98 C, followed by 2 cycles of: 20 seconds at 98 C, 30 seconds at Cell Reports Methods 1, 100070, September 27, 2021 e4 Report ll OPEN ACCESS 55 C and 40 seconds at 72 C, followed by x cycles of (x = number previously determined by qPCR): 20 seconds at 98 C, 30 seconds at 65 C and 40 seconds at 72 C, and terminated with 5 minute-incubation at 72 C. After PCR, the sample was purified by the addition of 0.7x AMPure XP magnetic beads and eluted in 10 ml of water. The quality of the final libraries was verified by ds-DNA concentration measurement using Qubit ds-DNA High Sensitivity assay (Thermo Fisher Scientific, USA) and BioAnalyzer (Agilent 2100, USA), following the manufacturer's instructions. If the quality of the samples was sufficient, the mRNA libraries were sequenced together with the protein libraries on with NextSeq500 (targeting 50 million reads per sample). The sequences of the primers used for the library preparations and sequencing are shown in Table S3.

Protein libraries preparation
After size separation of the protein and mRNA libraries, the protein library was amplified by IVT reaction. 20 mL of the product was mixed with 60 mL of IVT mix and amplified for 15 hours at 50 C. The IVT mix was composed of 1x reaction buffer, 10 mM each dNTP (ATP, CTP, GTP and UTP), and 10x diluted T7 RNA Polymerase mix (HiScribe TM T7 High Yield RNA Synthesis Kit, New England Biolabs, USA). After linear amplification, the sample was purified by the addition of 1.5x AMPure XP magnetic beads and eluted in 20 ml of water. 10 ml of the sample were stored at -80 C as post-IVT backup and the other 10 ml were mixed with 1 mL of 10 nM of each dNTP (Thermo Fisher Scientific, USA), 1 mL of 10 mM PE2-NNNN-Next1 primer (Biolegio, The Netherlands), and 1 mL of 10 mM PE2-NNNN-BioHash2 primer (Biolegio, The Netherlands) and reversely transcribed by incubation for 3 minutes at 70 C followed by cooling down the mixture on ice. Next, 1.5 mL nuclease-free water, 4 mL of 5x PrimeScript buffer (Takara, USA), 1 mL RNAse OUT (Thermo Fisher Scientific, USA) and 0.5 mL of 200 U/mL PrimeScript RT (Takara, USA) was added to the reaction mix that was incubated: at 30 C for 10 minutes, followed by 42 C for 1 hour and by 15 minutes at 70 C. After reverse transcription, the product was purified with 1.5x AMPure XP magnetic beads and eluted in 10 mL of nuclease-free water. Subsequently, the sample was PCR amplified after prior determination of the number of cycles needed by qPCR: 0.5 mL of purified material was mixed with 6.5 mL nuclease-free water, 10 mL of 2x Kapa HiFi HotStart PCR mix (Roche, Switzerland), 1 mL of 20X EvaGreen Dye (Thermo Fisher Scientific, USA), and 2 mL of PE1/PE2 primer mix (5 mM each) (Biolegio, The Netherlands). Then the PCR reaction was started by adding: 0.5 mL of nuclease-free water, 12.5 mL of 2X Kapa HiFi HotStart PCR mix (Roche, Switzerland), and 2.5 mL of PE1/PE2 primer mix (5 mM each) (Biolegio, The Netherlands) to the 9.5 mL of the sample. The PCR thermal cycling was as follows: 2 minutes at 98 C, followed by 2 cycles of: 20 seconds at 98 C, 30 seconds at 55 C and 40 seconds at 72 C, followed by x cycles of (x = number previously determined by qPCR): 20 seconds at 98 C, 30 seconds at 65 C and 40 seconds at 72 C, and terminated with 5 minute-incubation at 72 C. PCR amplified product was purified by addition of 1.2X AMPure XP magnetic beads, eluted in 10 ml water and its quality was verified by ds-DNA concentration measurement with Qubit ds-DNA High Sensitivity assay (Thermo Fisher Scientific, USA) and Bio-Analyzer (Agilent 2100, USA) characterization. If the quality of the samples was high enough, the protein libraries were sequenced together with the mRNA libraries on with NextSeq500 (targeting 50 million reads per sample). The sequences of the primers used for the library preparations and sequencing are shown in Table S3.
Bulk mRNA Sequencing data analysis RNAseq fastq files of samples were aligned to the human genome GRCh38 using HISAT2 v2.1.0. Number of reads was assigned to genes by using featureCounts v1.6.1. Reads mapped to genes were normalized using count per million method (CPM) implemented in edgeR package in R Bioconductor. To analyze genes that have the highest fold change upon different durations of stimulation with aIg (60, and 180 min), genes were ranked by the difference between their CPM normalized values and mean of CPM normalized values of this gene in the unstimulated control sample.
From FASTQ files to count tables After sequencing, sequence data from the NextSeq500 (Illumina) was demultiplexed. Two separate FASTQ files were generated for the mRNA and protein libraries. The quality of the sequencing data was evaluated using a FastQC tool (version 0.11.7, Babraham Bioinformatics). The Python script presented by Adrian Veres (https://github.com/indrops/indrops) was used with some modifications, (in particular, for the protein libraries analysis), to process the FASTQ files in order to generate an mRNA and a protein count table. Briefly, the FASTQ files were first filtered to verify that the reads had the correct structure, sufficient quality and complexity. In our setup, 70-80% of the reads were successfully passing the filtering process. Secondly, ''real'' cells were determined by identifying the cell barcodes having the highest number of reads and setting a threshold below which the barcodes with fewer reads are removed (i.e., cell encapsulated with two beads, a cluster of cells encapsulated with one bead, cells with degraded mRNA). Thirdly, reads were sorted according to their barcode of origin. Subsequently, the remaining reads were aligned using Bowtie with different settings for the mRNA and protein files. For the mRNA m = 10, n = 1, l = 15, e = 1000 as described by Adrian Veres and for the proteins m = 1, n = 1, l = 8, e = 75. With m = the maximum number of different alignments allowed per read, n = the number of mismatches allowed in the first l bases of the read and e = the maximum sum allowed of the quality values al all mismatched positions. The mRNA reads were aligned to the Homo Sapiens Genome FASTA file, while the protein reads were aligned to a self-written FASTA file in which were present the 8 bp Ab tag sequences. Finally, the reads were quantified to UMIFM (UMI-filtered mapped) counts and were aggregated in two separate count tables for the mRNA and for the proteins.