Advances in "-omics" technologies now provide scientists with the ability to probe human biology and biologic variance with high sensitivity at the single cell level. These approaches have been of particular benefit to studies of immune-mediated diseases in humans where deep profiling of peripheral blood and peripheral immune cells has provided insights to underlying pathobiology, unique biomarkers of disease, and biological variation. The ability to reliably discern the state of individual cells and discover true biologic heterogeneity in a complex system like peripheral blood requires that the effects of sample handling and storage not overwhelm those associated with the underlying biology. The logistical challenges related to collecting, processing, and shipping blood in human studies and particularly in clinical trials can make the application of these highly sensitive technologies a challenge. At present, there are no clear, published, multi-omic studies that provide practical guidance about how quickly peripheral blood samples need to be processed and which cells and cell pathways are most impacted by delayed processing.
Peripheral blood mononuclear cells (PBMC) are a workhorse of human immunology owing to the ease of collection and simplicity of cell isolation. Whole blood, from which PBMC are derived, is understood to be remarkably stable and previous work has established support for flexibility in sample handling for various whole blood assays, often as a way of managing the inherent challenges of human blood collection. Studies have shown that, when collected into anti-coagulant tubes to prevent clotting, whole blood left to "rest" was shown to be stable up to 24 hours at room temperature (Wu et al., 2017; Zini, 2014) and plasma profiling for common metabolites was similarly robust (van Eijsden et al., 2005; Zimmerman et al., 2012) .
Other data, however, demonstrated the potential for profound changes resulting from variability in sample handling, particularly in the context of untargeted assays such as those for transcriptomic analysis. Overnight delay of PBMC isolation from whole blood was shown to alter thousands of genes, in particular JUN, FOS, and the heat shock pathway (Baechler et al., 2004) , and even delays as short as four hours led to substantial changes, especially in immune-related gene expression (Barnes et al., 2010; Massoni-Badosa et al., 2020) . Particularly problematic are processing artifacts which impact a wide variety of genes and proteins related to immunity in health and disease, which can obscure the disease processes of interest (Dvinge et al., 2014) . Missing from previous analyses, however, are state-of-the-art assays such as single-cell RNA-sequencing that are becoming the standard for disease-based research and a comprehensive, multi-modal approach to understand the full complexity of ex vivo biology.
In an effort to clearly address these questions prior to initiating a series of multi-center clinical studies focused on human immunology, we performed deep, multi-omic profiling of human peripheral blood stored in anticoagulant for varying lengths of time before processing to plasma and peripheral blood mononuclear cells (PBMC) . This resource provides clear insights into the rapid changes related to delayed sample processing, elucidating the cells most altered, the cellular pathways most impacted, and the assays most affected. While flow cytometry did not reveal large-scale changes in cell type frequencies through an 18-hour delay, single-cell gene expression and high-plex plasma proteomics provide overwhelming evidence that cells of all types exhibit time-dependent activity that distorts the underlying biology. These changes are broad and dynamic, complicating the technical analysis of single-cell RNA-seq data and especially inferences of in vivo physiology from ex vivo assays. We propose the affected proteins and genes be carefully considered in any human biology study or clinical trial that uses blood and/or PBMCs to reflect in vivo biology, and that these findings may extend to blood-and immune-cell permeated tissue, as well. Finally, the baseline cytometry, proteomic, and transcriptomic data on 10 donors can serve as a highquality resource to accelerate human systems immunology research and provide the substrate to begin decoding these effects in existing and emerging studies.
To study the effects of delays in PBMC processing from whole blood we performed two similar but independent experiments ( Figure 1A) . In Experiment 1, we isolated PBMC or plasma from whole blood at 2, 4, 6, 8, and 18 hours after blood draw from healthy donors (n=3) or those diagnosed with systemic lupus erythematosus (SLE, n=3) . In Experiment 2, we assayed PBMC or plasma isolated from only healthy donors (n=4) starting at 2, 4, 6, 10, 14, and 18 hours after blood draw. In both experiments, the whole blood was held in the dark at room temperature prior to PBMC isolation by Ficoll gradient separation and PBMC were assayed after freeze/thaw by flow cytometry and 10x Genomics single-cell RNA-sequencing, and the plasma by Olink proteomics. Details of the samples, the assays used in each experiment, and any deviations are available in the Methods.
We started by assessing basic metrics from our technical processes. Cell yields from PBMC isolation, and their viability, did not exhibit time-dependent changes (Supplemental Figures  S1A-B) . Nor were the cells somehow more "fragile", as they exhibited similar post-thaw recovery and viability at any time of isolation (Supplemental Figures S1C-D) . This overall consistency was substantiated by flow cytometry profiling of the major immune cell types present in PBMC, most of which did not show noticeable changes throughout the processing delay time series ( Nicholson et al., 1984) , suggesting an increased recovery of low-density neutrophils by Ficoll gradient separation with longer delays in PBMC processing. In addition, more granular flow cytometry profiling in Experiment 2 showed decreases in plasmablasts, non-classical monocytes, and basophils ( Figures 1D-F) . Nevertheless, with the exception of rare neutrophils, time-dependent effects on the abundance of cell types were modest.
The minor effects seen by cytometry could result from subtle changes in PBMC isolation and PBMC survival, but little can be inferred of cellular activation or stress from the flow cytometry panels used in this study to quantitate cell types. To get a better understanding of cellular activity, we first assayed bulk PBMC from the six donors in Experiment 1 by targeted transcriptomics with the Nanostring nCounter platform. This clinical assay enumerates a panel of 594 pre-selected immune-and disease-relevant transcripts from a bulk cell sample. Despite well documented person-to-person variability in PBMC gene expression (Kaczorowski et al., 2017; Radich et al., 2004; Whitney et al., 2003) , Principal Component Analysis (PCA) indicated that PBMC processing delay accounted for a substantial part of the variability (Supplemental Figure S3) . These effects were reflected in unsupervised clustering of the data, where it was apparent that expression levels across nearly the entire targeted gene set are inverted at 18 hours post-blood draw compared to the previous time point (Figure 2A , box indicated by the solid line). A feature of this longitudinal data set is that gene regulation dynamics can be observed, revealing that in all six donors a subset of these genes underwent an earlier induction before returning to low levels (box indicated by the dashed line). Aggregate analysis highlighted common and progressive exaggeration of changes over the time course ( Figure 2B) , with 85 genes significantly induced (including IL20, DEFB4A, DUSP4, LILRB5) and 261 genes significantly down-regulated (including PD-1 (PDCD1), ST2 (IL1RL1), CX3CL1, and CCR2) at 18 hours relative to 2 hours post-draw ( Figure 2C -D and Supplemental Table S1 ). Without the context of the longitudinal data, these expression patterns could suggest physiologic biology that was not actually present in the host. Therefore, we conclude that while the cells themselves are generally robust to processing delays, their biological profile is greatly impacted.
Delayed processing results in substantial changes to the plasma proteome
The extensive changes recognized by bulk transcriptomics suggested that whole blood awaiting processing was an active and dynamic immune environment, and for numerous reasons we hypothesized there would be significant changes in plasma proteins as blood awaited processing: it has previously been reported that storage of whole blood led to increased levels of thrombospondin by eight hours (Kaisar et al., 2016) , immune cells are known to respond to hemostasis (Jenne et al., 2013) , temperature (Cho et al., 2010) and oxygen levels (Nizet and Johnson, 2009) , and many transcripts demonstrating increases in our bulk transcriptomics assay encode secreted proteins (Supplemental Table S1 ). We speculated that not only were cells likely undergoing an intrinsic response, but that the extrinsic environment was likely changing as well.
To test this hypothesis, we interrogated the 54 samples in our study for 1161 plasma proteins using a dual-antibody targeted assay (Söderberg et al., 2006) . We found nearly one-third of the 1161 proteins assayed showed relative abundance varying with delayed processing in a combined analysis of both experiments ( Figure 3A and Supplemental Table S2 ). While the number of proteins exhibiting significant change was relatively constant over the time course (350 proteins at 4 hours to 469 proteins at 18 hours post-draw), the magnitude of the changes exhibited a shift from 6 hours onward. Only 5 proteins were 1.5fold changed between 2 and 4 hours post-draw (LAT, CD40L, EGF, PDGFβ, and SDC4), while by 6 hours 27 proteins had changed by more than 1.5-fold and by 18 hours 69 proteins had changed by at least 1.5-fold, relative to the 2-hour time point. This suggested an increasing activity within the samples as they awaited processing.
Some of the most highly affected protein levels suggest complex and coordinated activity. For example, AnxA1 and S100A11 directly interact on early endosomes, and EGFR can both lead to phosphorylation of and be degraded by AnxA1 (Poeter et al., 2013) . ANXA1, S100A11, and EGF protein are all increased in the plasma over time. More striking was that the number of platelet-related proteins significantly changed over time. In agreement with previous work (Khan et al., 2006) , one of the earliest and most increased proteins was CD40LG (Figures 3B-D and Supplemental Table S2 ), a product of platelet activation (Henn et al., 1998) . PDGFβ, ANGPT1, STK4, STX8, OSM, and SDC4, were also rapidly and significantly increased. These proteins are all known to be involved in platelet biology (Andrae et al., 2008; Beck et al., 2017; Golebiewska et al., 2015; Londin et al., 2014; Tanaka et al., 2003) , perhaps indicative of regulation or dysregulation of platelet activation in the waiting blood despite the presence of anti-coagulant. Fewer proteins exhibited decreases in abundance, with APLP1 the most striking. APLP1 can be processed by gamma-secretase but, unlike the homologs APP and APLP2, it was suggested that cleavage of APLP1 could occur without ectodomain shedding (Schauenburg et al., 2018) , potentially sequestering the protein intracellularly. MMP7 also exhibited a significant decrease over the time course and was shown to be a target of platelet activation (Yang et al., 2020) and platelet-derived CXCL4 (Erbel et al., 2015) . These data demonstrate that many of the strongest changes in plasma proteins in whole blood ex vivo relate to platelet biology.
While the majority of cell-type frequencies were essentially stable throughout the blood processing delay (Figure 1) , plasma proteomics showed the signals available to those cells were dynamic and the corresponding bulk transcriptomics indicated similar changes in gene regulation, as might be expected from an immune system evolved to respond to extracellular 7 signals. To better understand the impact of blood processing delay on distinct cell types we performed droplet-based single-cell RNA-sequencing on both experiments. In total, we collected and analyzed 450,189 high quality singlet cells with an average median features per cell of 1942 (Supplemental Figure S4 and Supplemental Table S3 ), for a total of more than 800 million unique RNA molecules quantitated by single-cell RNA sequencing.
We assessed the overall effect of delayed blood processing across all cell types by aggregating all samples in each experiment separately and summarizing the highdimensional data set in two dimensions using tSNE (Figures 4) . With this visualization, it was clear that by 18 hours post-draw (purple) there was a profound shift in the global gene expression pattern per cell across multiple clusters, with the time effect dominating interdonor variability. While this effect was clear and pervasive at 18 hours in Experiment 1, we were concerned that the irregular time spacing (PBMC isolation at 2, 4, 6, 8, and 18 hours post-draw) was a confounder in the way the data clustered. We therefore changed the experimental design of Experiment 2 (PBMC isolation at 2, 4, 6, 10, 14, and 18 hours postdraw) so the later time points occurred at more regular intervals. The time-dependent divergence from the native transcriptional state was even more pronounced in Experiment 2, with a progressive distancing apparent at 10 hours (light blue), 14 hours (dark blue), and 18 hours (purple) post-draw. These effects were not restricted to a limited number of cell types, such as monocytes and dendritic cells (DC), but were apparent in all labeled cell types, even those thought to be quiescent, such as naive T cells (Sprent and Surh, 2011) . These data demonstrate that analysis of global gene expression profiles in human PBMC can be confounded by pre-analytical process artifacts.
We were interested to more deeply understand the effect of blood processing delay on individual cell types to better define the changes on functionally distinct immune cell populations. To highlight cell type specific changes, we isolated the tSNE coordinates of the clusters mapping to each cell type (Figures 5A) and quantified the occupancy of each cluster over the time course. In this approach, we infer that each cluster represents a distinct transcriptional state of a given cell type, such that a change in cluster composition within a cell type indicates changing transcriptional states. For example, most CD14 + monocytes start from a single, dense cluster at the earliest time point which disintegrates concurrently with the appearance of three new clusters by 18 hours, with the major transition occurring between 4 and 6 hours post-blood draw ( Figure 5B) . A similar pattern was observed for CD4 + memory T cells, CD8 + memory T cells, and double-negative T cells. In contrast, the originating clusters of CD8 + naive T cells and CD56 low NK cells persisted throughout the time course despite the appearance of new clusters, and CD4 + naive T cells exhibited some clusters persisting, some disappearing, and some appearing throughout the time course. Though the effects on different cell types resulted in different patterns of changing transcriptional profiles, for nearly all cell types changes were most apparent between the 4hour and 6-hour time points ( Figure 5C ). This was true of cell types thought to be most responsive to environmental signals (e.g. monocytes) and for those thought to be quiescent (e.g. CD4 + naive T). From these data we conclude that measurements obtained from PBMC processed more than four hours after blood draw are increasingly at risk for pre-analytical artifacts.
Reflecting the changes in the global gene expression profile observed through unsupervised clustering, quantification of differentially expressed genes (DEG) relative to the 2-hour time point demonstrated a trend of an increasing number of DEG across cell types (Supplemental Figure S5 and Supplemental Table S4 ). Once again, the number of changes was muted when PBMC processing started at 4 hours but emerged following the 6hour time point. As expected, the specific transcriptional changes exhibiting time-dependent regulation were varied across cell types, but generally reflected an elevated activation status and included many shared transcripts, such as NFκB pathway genes (NFKBIA, REL, TNFAIP3), CD83, JUN family members, HIF1A, and EIF1. But coordinated regulation was also apparent, such as induction of MAP3K8, IRF2BP2, TLE3, and MIR22HG in both monocyte subsets and classical dendritic cells, and SBDS and RBM38, both RNA binding proteins, in all eight lymphocyte subsets.
Down-regulated genes exhibited similar trends. DDX17, RIPOR2, TRAF3IP3, UCP2, and CD53 were diminished across many of the 14 cell types classified in the scRNA-seq data, while DYNLL1, IFI16, POU2F2, CALM2, FYB1, and FGL2 were coordinately down-regulated in monocytes and dendritic cells in both experiments (Supplemental Table S4 ). We were particularly interested in the GIMAP (GTPase of the immunity-associated protein) family, a collection of eight sequence-related genes implicated in pro-and anti-apoptotic functions, primarily through studies in mouse T cells (Filén and Lahesmaa, 2010) . GIMAP4 and GIMAP7 were strongly down-regulated in nearly all cell types, with the exception of B cells and plasmacytoid dendritic cells, and GIMAP1, GIMAP6, and GIMAP8 also exhibited decreasing expression in various cell types over the time course. In fact, GIMAP4, GIMAP7, and GIMAP1 were among the genes that had the greatest drop in cellular percentage detection, irrespective of expression level. These data suggest survival or apoptosis may be key processes impacted by delayed PBMC processing.
To more comprehensively assess the changes occurring through the processing delay we performed pathway analysis by comparing the dominant cluster late in the time series and the dominant cluster early in the time series within a given cell type. For this analysis we selected CD14 + monocytes and CD4 + memory T cells as examples of myeloid and lymphoid cell types frequently interrogated in systems biology studies. Through this analysis, 18 of the top 22 pathways in CD14 + monocytes and 23 of the top 26 pathways in CD4 + memory T cells enriched in the late clusters included NFκB and/or JUN family members (Supplemental Table S5 ), indicating the emergence of an activated phenotype within these cell types. Also of note was the elaboration of pathways in both cell types related to HMGB1, hypoxia, Th17 activation, the NFκB pathway, IL-6 signaling, and apoptosis. Other pathways of note include two AHR pathways and TREM1 signaling induced in CD14 + monocytes, and the mTOR and unfolded protein response in CD4 + memory T cells. These pathways variously relate to survival and apoptosis, inflammation, and sensing of extracellular signals that potentially confound the native ex vivo biology. Though we cannot determine whether this shift in cluster occupancy and gene expression program represents a direct precursor-progeny relationship or whether it results from the outgrowth of a distinct rare subset of cells within each population, in either case it is clear that cells labeled as CD14 + monocytes or CD4 + memory T cells have undergone a time-dependent shift in their transcriptional profile. On the whole, the single-cell transcriptomics data demonstrate that cells are not "resting" in whole blood in the time between blood draw and PBMC processing, but rather are subject to and participating in a complex and active immune environment.
A unique advantage of these data is the ability to use them in cross-modal analyses. All data sets from a given donor were acquired from the same blood draw, eliminating subclinical immunity, circadian, diet, sleep, seasonal allergy, and a host of other environmental variables as confounders for intra-donor longitudinal analysis. In addition, targeted bulk transcriptomics and single-cell transcriptomics were performed from the same PBMC isolation per time point. Therefore, these single-cell RNA-seq data are ideal for deconvoluting the bulk transcriptomics of PBMC. One striking example is the T cell costimulatory molecule ICOS ligand (encoded by ICOSLG). Targeted bulk transcriptomics indicates a modest log2 fold change of 1.7-2.3 across the six donors in Experiment 1 ( Figure 6A ). Single-cell RNA sequencing, however, identifies induction of ICOSLG not only in myeloid cells such as conventional dendritic cells and plasmacytoid dendritic cells but also in B cells, with the highest expression in naive B cells ( Figure 6B ). Another interesting example is the secreted plasminogen activator urokinase, encoded by the gene PLAU, which was identified in the bulk transcriptomics data for a moderate induction by the 18-hour time point of Experiment 1. Single-cell analysis indicates the majority of this induction occurs in plasmacytoid dendritic cells, with little expression in other cell types. Urokinase is localized to the extracellular plasma membrane by association with PLAUR (Ellis and Danø, 1991) , which was likewise induced in bulk transcriptomics and shown to be expressed by monocytes and dendritic cells in our data, thus establishing a potential source and reservoir for potentiation of thrombolytic activity. These examples highlight the potential to explicitly decode bulk transcriptomic data sets of PBMC using our six matched longitudinal data sets.
More intriguing is the potential for synergy with cross-modal data. Whole blood left on the benchtop to be processed is a closed system and observations of concordant changes in plasma protein and the corresponding transcript or cognate receptor transcript could be used to build testable hypotheses. To find such protein-transcript pairings, we cross-referenced the most varied plasma proteins with the corresponding scRNA-seq data, highlighting a number of interesting targets. One, the TNF superfamily member LIGHT (encoded by TNFSF14), is a membrane and secreted ligand of the receptors HVEM, LTβR, and the soluble decoy DcR3. LIGHT potentiates T cell responses (Shaikh et al., 2001) , has been reported in the setting of autoimmunity (Herro et al., 2015; Kotani et al., 2012) , and is actively pursued as an anti-cancer therapy (Skeate et al., 2020) . In our proteomics data, detectable LIGHT protein in the plasma fraction increased over the blood processing delay time course of both experiments ( Figure 7A ), resulting in a mean 2.49 fold increase at 18 hours post blood draw relative to 2 hours, the 20 th most-induced protein in that comparison.
Simultaneously, single-cell RNA sequencing showed the highest expression of LIGHT in CD8 + memory α β
α β T and many γ δ T cells) ( Figure 7C ). One of the receptors of LIGHT, HVEM (encoded by TNFRSF14), was broadly expressed but the other major receptor LTBR was restricted primarily to the myeloid compartment. At the same time, the decoy receptor DcR3 (encoded by TNFRSF6B), capable of negatively regulating LIGHT (Wroblewski et al., 2003) , was expressed mostly in lymphocytes, especially by double-negative T cells, which showed a trend toward increasing transcript levels concordant with increased levels of LIGHT plasma protein.
Another intriguing participant in the complex ex vivo blood environment is BACH1. BACH1 is a BTB/POZ transcription factor that heterodimerizes with small MAF proteins, leading to the regulation of oxidative stress gene targets bearing MAF recognition elements (Zhang et al., 2018) . One such target, HMOX1, can be rapidly induced in monocytes by inflammatory stimuli (Yachie et al., 2003) , and in fact we recognized induction and expression of HMOX1 in CD14 + and CD16 + monocytes over the processing delay in our data ( Figure 7D) . BACH1 was identified in our data sets by an increase in plasma protein levels ( Figure 7B ) and has previously been recognized as a regulator of HMOX1 in human monocytes (Miyazaki et al., 2010) . Consistent with the literature, we found the highest levels of BACH1 transcript were expressed by monocytes, especially early after blood draw ( Figure 7D) . However, whereas nearly all other populations exhibited little to no expression of BACH1 throughout the time course, both naive and non-naive B cells demonstrated a sharp induction of BACH1 gene expression. A role for BACH1 in B cells has not been reported aside from potential redundancy with BACH2 during B cell development (Itoh-Nakadai et al., 2014) . Thus, this observation represents a potential novel aspect of B cell biology under conditions of cellular stress and highlights the potential for these aligned data sets not only to reveal putative sources of activation-induced proteins, but to suggest possible cellular networks of communication.
1 0
The human immune system can be biopsied through the collection of peripheral blood samples, making blood one of the most amenable human tissues for research studies. As a 'sense and respond' organ, the immune system tailors its activity to a wide variety of subtle environmental cues, including ex vivo manipulation (Kelley et al., 1987) . We used a longitudinal multi-omics approach to better understand artifacts from process variability in the study of whole blood components, providing three important contributions to the study of the human immune system. First, we showed that a delay in processing venipuncture blood samples has profound consequences on immune cells. Global effects on the native biology were identified in multiple assays, in all 10 donors of the study, and in all cell types we identified. Second, our multi-modal and longitudinal data sets constitute over 10 billion measurements of immune activity in a closed system. These data can be mined for correlative gene-gene, gene-protein, ligand-receptor, and pathway interactions from which testable hypotheses can be derived. Furthermore, as these activities relate to potentially confounding pre-analytical artifacts, these data may enable a qualification of results from other studies that are subject to variable or sub-optimal processing conditions. Third, the 2hour time points provide a large, high-quality, multi-modal data set of 10 human donors. These 2-hour data are temporally close to the in vivo biology, consisting of flow cytometry for major PBMC populations, single-cell RNA-sequencing of ~8000 quality cells per donor, targeted plasma proteomics of >1000 targets, and targeted bulk transcript analysis of >600 targets on six donors. Thus, this resource provides various opportunities to investigate a multi-modal dataset derived from the components of freshly isolated and ex vivo aged human blood.
While the data from typical flow cytometry implied that commonly used cell surface proteins and the cells themselves are stable, both the targeted as well as whole transcriptome studies confirm large-scale changes in gene regulation. After an 18-hour delay in blood processing, transcriptional states in all cell types were fundamentally altered. Given the length of time, this might be expected. But more concerning are the changes occurring after just 4 and 6 hours post-draw, as evidenced by the changes in cellular state apparent through single-cell RNA-sequencing. Our findings implicate an up-regulation in inflammation and stress-activation (NFKBIA, HIF1A, MAP3K8), perhaps in response to deteriorating environmental conditions (Paardekooper et al., 2018; Rius et al., 2008) . Under these conditions, the appearance of new clusters in tSNE visualization of gene expression data likely results from a large number of cells altering their gene expression profiles ex vivo rather than outgrowth of more rare pre-existing cells faithfully representing an in vivo activation state. If true, this fact underscores the challenges of using tSNE/UMAP classification of the data as an analysis mode, as cluster occupancy can be confounded by blood processing delay. The preferred alternative is to identify cellular identity on a per-cell basis first, and then compare the transcriptomic state of equivalent cell types based on the experimental variable (e.g. disease state); in this case the sample processing metadata, such as blood processing delay, must be carefully tracked and corrected for in the analytical model.
As predicted from our gene expression data, the plasma proteome exhibited dramatic changes in concert with the PBMC transcriptomics data. Efforts to understand the effect of processing delay on various plasma and serum analytes have existed for some time (Ono et al., 1981) , though typically limited to individual proteins or metabolites of known clinical interest. We have generated a unique public human proteomic data set, consisting of over 62,000 data points linked to PBMC flow cytometry and transcriptomics. The clearest signal in the targeted proteomics relates to platelet biology. This includes CD40L, PDGFβ, ANGPT1, STK4, and others, supportive of a hypothesis of early and ongoing platelet activation. Blood for the proteomics studies was drawn into K2-EDTA tubes, the recommended anticoagulant for the assay and consistent with common practice. EDTA has previously been reported to 1 1 elicit aggregation of platelets in vitro (Pegels et al., 1982) . Of note, platelet aggregation in whole plasma is also not inhibited by heparin (Saba et al., 1984) , such that the platelet response to blood draw or hemostasis may elicit or potentiate divergence from native biology. Therefore, our data suggests that some of the most severe changes in plasma proteomics result from the basic process of obtaining the sample. This thrombotic-like activity has the potential for cascading effects on PBMC transcriptional activity and may be a key source of in vitro artifacts. However, it is clear that such changes are exacerbated by time delay, thus highlighting the need for rapid and consistent processing.
In order to mitigate these confounding effects, the optimal solution is to process the samples consistently and with as little delay as possible. Previous studies on bulk PBMC have suggested that storing samples on ice may mitigate the effects of delayed blood processing (Goods et al., 2018) , and so to test this hypothesis we kept some of the whole blood samples at 4°C in a parallel experiment (data not shown). We found this mitigation effort to have mixed results, with an unexpected decrease in PBMC yield and hemolysis occurring in some samples. Therefore, our data indicates processing within 4 hours should become standard practice for assays using human blood components, particularly those employing deep transcriptional or proteomic profiling. This may not be possible in all cases. In these situations, our data provide a comprehensive set of landmarks to qualify the analysis for processing delay artifacts.
Embedded within the larger resource provided here, the 2-hour cytometry, proteomic, and transcriptomic data sets from 10 donors is a substantial resource in its own right, with minimal and controlled processing delay. The available public cytometry and single-cell RNA-seq data on human PBMC continue to grow in size and complexity, but to our knowledge no study integrates these two modalities with deep proteomics at this scale. Furthermore, though blood processing delay results in immune activity that may misrepresent or obscure the true biology being studied, we suspect that in many cases this 'artifactual' biology echoes physiologic processes. This potential was highlighted by the broad down-regulation of multiple GIMAP family genes in many cell types. As these are known to relate to apoptosis in T cells (Filén and Lahesmaa, 2010) , this observation is interesting. However, our data also opens the possibility of a similar role in myeloid cells, for which there is little existing data (Hellquist et al., 2007; Krücken et al., 1997) . Therefore, the multi-modal data set is a substrate for generating testable hypotheses, so long as the limitations of an ex vivo artificial system are fully accepted. In particular, blood flow has ceased and cells may aggregate by settling, there is no regulation of blood gases, and temperature is well below physiologic range. When approached with an understanding of the experimental conditions employed, we believe these data can be used productively, for example in the study of transcriptional co-regulation, correlations among proteins and between proteins and genes, as well as novel associations of genes with cell types.
Recently, Massoni-Badosa et al. published a complementary effort on the effect of blood processing delay (Massoni-Badosa et al., 2020), with differences in the assays, cell numbers, and analyses performed. Our overall conclusions align with theirs, in the clear risks to data interpretation from processing delay and non-uniformity, though we did not recognize a global or cell type-specific reduction in transcript levels (Supplemental Figure S4) . Nonetheless, the practical realities of obtaining clinical samples will ultimately dictate how and when they are processed, and the data generated should be analyzed with an appreciation for the unknown and/or variable amount of delay in processing. While this naturally applies to myriad ongoing and historical clinical studies, it is also keenly relevant to public health imperatives such as the SARS-CoV2 pandemic, for which process standardization may be outweighed by the benefit of sample acquisition. Accounting for such practicalities in a data-driven manner will provide broad and on-going 1 2 benefit for improving the caliber of putative therapeutic targets, by lowering confidence in those subject to ex vivo regulation and thereby increasing confidence in others. We hope these data will also spur a renewed focus on process standardization and optimization to ultimately improve data quality and physiologic relevance at the point of data generation. With over 1 billion data points across several -omic platforms, our resource offers a rich data repository upon which to qualify existing data, benchmark and train subsequent studies, and derive novel hypotheses. 
The authors have nothing to declare. Figure S3 and Supplemental Table S1 . (B-C) All proteins in Experiment 1 (B) and Experiment 2 (C) were plotted for fold change and adjusted p-value (Benjamini-Hochberg FDR) as calculated using a GLMEM model, and colored by time point.
(D) A selection of proteins showing the greatest fold change in combined analysis were plotted by time point and shaded regions indicate the 95% confidence interval. Where both experiments identified a significant fold change for a given protein, but those fold changes differed between experiments according to the GLMEM model, two lines are shown (dashed for Experiment 1, dotted for Experiment 2); otherwise a single line (solid) representing data from both experiments is shown. Table S2 .
The single-cell RNA-sequencing normalized gene expression matrices from each sample were identified for donor and time point and aggregated per experiment (Experiment 1: donors 2-7; Experiment 2: donors A-D). The multidimensional data was displayed in two dimensions using tSNE, colored as indicated. The data are overlaid in sequence with the latest time point (18hr) on top, obscuring some data points. Additional details can be found in the Methods.
Individual cells were assigned a cell type label independently of clustering or tSNE visualization.
(C-F) tSNE visualizations colored by time point for Experiment 1 (C) and Experiment 2 (D). The tSNE maps in (C) and (D) were split out by donor for Experiment 1 (E) and Experiment 2 (F). Figure S4 and Supplemental Table S3 .
Single-cell RNA-sequencing data from Experiment 2 was analyzed as in Figure 4 (C) The number of cells occupying each cluster were enumerated and the counts (y-axis) plotted per cell type relative to the time after blood draw (x-axis) for each donor. Colors match the Louvain cluster colors in (B). Thicker lines represent the mean for all donors and thinner lines represent single donors. The 2-4 hour range highlighted in pink shows a period of apparent stability prior to more extreme divergence in cluster occupancy after 4 hours.
See also Supplemental Figure S5 and Supplemental Tables 4, 5, and 7.
(A) Transcripts from bulk transcriptomics data analyzed as in Figure 2 were selected based upon increasing expression over time and consistent dynamics across all six donors in Experiment 1. Charts show normalized counts for hours since blood draw for the individual donors, indicated by color.
(B) Single-cell RNA-seq data (analyzed as in Figure 4 ) from all 10 donors in both experiments were queried for transcripts selected in (A) and plotted by cell type. Charts show the normalized and scaled RNA count (y-axis) for hours since blood draw (x-axis). Individual donors are indicated by color.
(A-B) Plasma proteins from proteomics data analyzed as in Figure 3 were selected based upon increasing expression over time and across all 10 donors in both experiments. Normalized protein expression (y-axis) over time since blood draw (x-axis) is shown, with individual donors indicated by color.
(C-D) Single-cell RNA-seq data (analyzed as in Figure 4 ) from all 10 donors in both experiments were queried for transcripts encoding TNFSF14 and its receptors (C) or BACH1 and a selection of its transcriptional targets related to oxidative stress (D), and plotted by cell type. Charts show the normalized and scaled RNA count (y-axis) for hours since blood draw (x-axis). Individual donors are indicated by color.
Supplemental Table S1 (C-D) 5x10 6 PBMC isolated as in panels A-B and frozen in liquid nitrogen were thawed and assessed for yield (C) and viability (D) prior to downstream assays. "% recovery" was calculated as the recovery of live PBMC divided by 5x10 6 .
(E-F) Flow cytometry gating schemes for Experiment 1 (E) and Experiment 2 (F). Orange labels indicate gates used to determine population frequencies. The gate labeled "cleanup" in Experiment 2 was used to remove dye aggregates. Table S1 -Nanostring bulk transcriptomics principle components analysis PBMC in Experiment 1 were prepared from whole blood at various times after blood draw and assayed by targeted transcriptomics in bulk using Nanostring nCounter. The normalized counts were analyzed by principle components analysis and colored by donor (A) or by time point (B).
Single-cell RNA-sequencing technical metrics were compiled for Experiment 1 (A-D) and Experiment 2 (E-H) and categorized by cell type, as assigned by Seurat-based reference alignment (see Methods for details).
The number of quality singlets (A, E), mean mitochondrial DNA (B, F), features (C, G), and UMI (D, H) were quantified per cell type and colored by donor. 
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Adam Savage (adam.savage@alleninstitute.org).
No materials were generated for this study.
Flow cytometry .fcs data files are available on request. Nanostring bulk transcriptomics and Olink proteomics data are provided in Supplemental Tables S1 and S2, respectively. All single-cell RNA-sequencing data have been deposited at GEO (GSE156989).
Blood samples were obtained from healthy (no diagnosis of disease) and non-matched donors with systemic lupus erythematosus (SLE) (Supplemental Table S6 ), from the Benaroya Research Institute (Seattle, WA) or Bloodworks Northwest (Seattle, WA) through protocols approved by the relevant institutional review boards. Blood was drawn into BD NaHeparin vacutainer tubes (for PBMC; BD #367874) or K2EDTA vacutainer tubes (for plasma; BD #367863) and, upon arrival at the processing lab, all NaHeparin tubes for each donor were pooled into a sterile plastic receptacle to establish one common pool for all time points. PBMC isolation and plasma processing was started at 2, 4, 6, 8, and 18 hours postdraw (Experiment 1) or 2, 4, 6, 10, 14, and 18 hours post-draw (Experiment 2). Due to a PBMC processing error, flow cytometry and single-cell RNA-sequencing was not available for the 6-hour time point from donors A and B. This resulted in 30 samples of bulk transcriptomics (Experiment 1 only), 52 samples of flow cytometry and single-cell RNAsequencing, and 54 samples of plasma proteomics data.
For PBMC isolation, at each time point the pool of blood was gently swirled until fully mixed, about 30 times, and a volume of blood was removed and combined with an equivalent volume of room temperature PBS (ThermoFisher #14190235). PBMC were isolated using one or more Leucosep tubes (Greiner Bio-One #227290) loaded with 15 mL of Ficoll Premium (GE Healthcare #17-5442-03) to which a 3 mL cushion of PBS had been slowly added on top of the Leucosep barrier. The 24-30 mL diluted whole blood was slowly added to the tube and spun at 1000xg for 10 minutes at 20°C with no brake. PBMC were recovered from the Leucosep tube by quickly pouring all volume above the barrier into a sterile 50 mL conical tube. 15 mL cold PBS+0.2% BSA (Sigma #A9576; "PBS+BSA") was added and the cells were pelleted at 400xg for 5-10 minutes at 4-10°C. The supernatant was quickly decanted, the pellet dispersed by flicking the tube, and the cells washed with 25-50 mL cold PBS+BSA. Cell pellets were combined, if applicable, the cells were pelleted as before, supernatant quickly decanted, and residual volume was carefully aspirated. The PBMC were resuspended in 1 mL cold PBS+BSA per 15 mL whole blood processed and counted with a ViCell (Beckman Coulter) using VersaLyse reagent (Beckman Coulter #A09777) or with a Cellometer Spectrum (Nexcelom) using Acridine Orange/Propidium Iodide solution. PBMC were cryopreserved in Cryostor10 (StemCell Technologies #07930) or 90% FBS (ThermoFisher #10438026) / 10% DMSO (Fisher Scientific #D12345) at 5x10 6 cells/mL by slow freezing in a Coolcell LX (VWR #75779-720) overnight in a -80°C freezer followed by transfer to liquid nitrogen.
For plasma isolation, the K2-EDTA source tube was gently inverted 10 times and the appropriate volume of whole blood was extracted using an 18-gauge needle and syringe, and transferred to a similar plastic tube with no additives (Greiner Bio-One #456085). The blood was centrifuged at 2000xg for 15 minutes at 20°C with a brake of 1, and 80-90% of the plasma supernatant was removed by careful pipetting for immediate freezing at -80°C. Plasma was assayed after the first freeze/thaw.
PBMCs were removed from liquid nitrogen storage and immediately thawed in a 37°C water bath. Cells were diluted dropwise with 37°C AIM V media up to a final volume of 10 mL. A single wash was performed in 10 mL of PBS+BSA. PBMCs were resuspended 2 mL in PBS+BSA, counted on a ViCell or Cellometer Spectrum, as above, and 1x10 6 cells were incubated sequentially or together with Human Trustain FcX (BioLegend #422302) and Fixable Viability Stain 510 (BD #564406), on ice according to manufacturer's instructions. Cells were washed with PBS+BSA and stained with a master mix cocktail of antibodies on ice for 25-30 minutes. Experiment 1 and Experiment 2 were assayed by flow cytometry using a 12-target panel and a 24-target panel, respectively (Methods/Antibodies). Cells were washed with PBS+BSA and fixed with 4% paraformaldehyde (Electron Microscopy Sciences #15713) for 12-15 minutes at room temperature. Cells were washed with PBS+BSA, resuspended in PBS+BSA, and collected on a BD Symphony cytometer. After compensation, cytometry data was pre-processed to remove unrepresentative events due to instrument fluidics variability (time gating), to exclude doublets (by FSC-H and FSC-W), and to exclude cells exhibiting membrane permeability (live/dead gating) prior to quantification using BD FlowJo software. Gating examples are provided in Supplemental Figure 1E -F. The percent change in frequency was calculated per cell type to quantify the observed cell frequency changes along the variable of time, along with the corresponding change in median frequencies across all samples.
For Experiment 1, 2x10 5 fresh PBMC were pelleted, resuspended in Qiagen buffer RLT (Qiagen #79216), and quick frozen immediately at -80°C for assay by the Nanostring nCounter platform, performed as a fee-for-service by Nanostring using their standard protocols for the nCounter Gene Expression -Hs Immunology v2 CodeSet assay.
Plasma samples were assayed using the Olink Proximity Extension assay, run on the Fluidigm Biomark system. They were submitted to Olink Boston (Experiment 1) or Olink Sweden (Experiment 2). Patient samples were distributed evenly across plates, and all timepoints per patient were run on the same plate, with randomized well locations. Samples were assayed using the Olink Discovery Assay which encompasses all proteins across 13 panels ( 3311) ). Data were normalized by Olink using an Olink-provided plasma plate bridging control, three positive controls, and three background controls. Samples reported below the stated Limit of Detection were removed from our analysis.
PBMCs were thawed as described in the Methods for flow cytometry and the same vial was used for flow cytometry and single-cell RNA-sequencing. Single-cell RNA-seq libraries were 2 1 generated using the 10x Genomics Chromium 3' Single Cell Gene Expression assay (#100075 or #1000121) and Chromium Controller Instrument according to the manufacturer's published protocol. 16,000 cells from each PBMC sample were loaded into a separate Chromium Single Cell Chip B (10x Genomics #1000073) well targeting a recovery of 10,000 cells. Gel Beads-in-emulsion (GEMs) were then generated using the 10x Chromium Controller. The resulting GEM generation products were then transferred to strip tubes and reverse transcribed on a C1000 Touch Thermal Cycler programmed at 53°C for 45 minutes, 85°C for 5 minutes, and a hold at 4°C. Following the reverse transcription incubation, GEMs were broken and the pooled single-stranded cDNA fractions were recovered using Silane magnetic beads (Dynabeads MyOne SILANE #37002D). Purified barcoded, full-length cDNA was then amplified with a C1000 Touch Thermal Cycler programmed at 98°C for 3 minutes, 11 cycles of (98°C for 15 seconds, 63°C for 20 seconds, 72°C for 1 minute), 72°C for 1 minute, and a hold at 4°C. Amplified cDNA was purified using SPRIselect magnetic beads (Beckman Coulter #22667) and a 1:10 dilution of the resulting cDNA was run on a Bioanalyzer High Sensitivity DNA chip (Agilent Technologies #5067-4626) to assess cDNA quality and yield. A quarter of the cDNA sample (10 ul) was used as input for library preparation. Amplified cDNA was fragmented, end-repaired, and A-tailed is a single incubation protocol on a C1000 Touch Thermal Cycler programmed at 4°C start, 32°C for minutes, 65°C for 30 minutes, and a 4°C hold. Fragmented and A-tailed cDNA was purified by performing a dual-sided size-selection using SPRIselect magnetic beads (Beckman Coulter #22667). A partial TruSeq Read 2 primer sequence was then ligated to the fragmented and A-tailed end of cDNA molecules via an incubation of 20°C for 15 minutes on a C1000 Touch Thermal Cycler. The ligation reaction was then cleaned using SPRIselect magnetic beads (Beckman Coulter #22667). PCR was then performed to amplify the library and add the P5 and indexed P7 ends (10x Genomics #1000084). PCR was performed on a C1000 Touch Thermal Cycler programmed at 98°C for 45 seconds, 13 cycles of (98°C for 20 seconds, 54°C for 30 seconds, 72°C for 20 seconds), 72°C for 1 minute, and a hold at 4°C. PCR products were purified by performing a dual-sided sizeselection using SPRIselect magnetic beads (Beckman Coulter #22667) to produce final, sequencing-ready libraries. Final libraries were quantified using Picogreen and their quality was assessed via capillary electrophoresis using the Agilent Fragment Analyzer HS DNA fragment kit and/or Agilent Bioanalyzer High Sensitivity chips. Libraries were sequenced on the Illumina NovaSeq platform using S4 flow cells. Read lengths were 28bp read1, 8bp i7 index read, 91bp read2.
The Nanostring platfrom was run on all timepoints from Experiment 1 (Subjects 2-7, Timepoints 2,4,6,8,18 hours). The data were normalized to Nanostring-provided control genes using the NanoStringNorm package (Waggott et al., 2012) and log-transformed (base 2). Genes with zero expression in more than 20% samples were removed from further analysis. Generalized linear mixed effect models (GLMEM) were fit to the expression data of individual genes, controlling for the fixed effects of gender (male vs. female) and binary processing timepoints (4, 6, 8, 18 hours vs. 2 hour-baseline) and subject-level random effects to account for individuals' inherent variability, using Gaussian family distribution. All the models were built in the R library glmmTWB version 1.0.1. The log2 (fold change from the baseline) was evaluated by the corresponding time-related effect estimates. Multiple tests correction was applied to the p-values to control the false discovery rate (FDR) using Benjamini & Hochberg procedure (Benjamini and Hochberg, 1995) .
Plasma from Experiment 1 and Experiment 2 were submitted to Olink (Uppsala, Sweden) in two separate batches respectively. All 2-hour samples from Experiment 1 were included in both batches as bridging controls. Proteins below the reported Limit of Detection or that exhibited more than 20% missingness or ranked in the lowest 20% among all proteins in expression (as defined by the sum of the corresponding normalized expression values) were disregarded from further evaluation. As recommended by Olink, expression values for each protein in Experiment 2 were subtracted by the corresponding protein-specific normalization factor, which was the median of the corresponding pair-wise expression difference as measured on the bridging controls between the two experiments.
A protein-specific GLMEM was fit to each normalized protein, controlling for the fixed effects of gender (male vs. female), age (35-55 years old vs. 25-35 years old), disease status (SLE vs. healthy), study id (Experiment 2 vs. Experiment 1) and continuous processing time (hour) and accounting for an individual's inherent variability through subject-level random effects. To capture the study specific non-linear temporal effects, the cubic, quadratic or linear effect of processing time and its interaction with study id were evaluated in a step-down approach using a likelihood ratio test. The dispersion was modeled as a function of plate, disease status and study id to allow non-constant variability. All the models were built in the R library glmmTWB version 1.0.1. We reported the average log2 fold-change in proteomic expression between later processing time points and the first time point (2 hour) along with the corresponding two-sided 95% CIs. Wald tests were conducted to assess significance of the change estimates and associated p-values. When a significant difference was found in the change estimates between the two studies due to experimental conditions unaccounted for, we reported study-specific change estimates. Multiple tests correction was applied to the pvalues to control the FDR using Benjamini & Hochberg procedure (Benjamini and Hochberg, 1995) .
Binary Base Call (BCL) files were demultiplexed using the mkfastq function in the 10x Cell Ranger software (version 3.1.0), producing fastq files. Fastq files were then checked for quality (FastQC version 0.11.3) and run through the 10x Cell Ranger alignment function (cell ranger count) against the human reference annotation (Ensembl GRCh38). Mapping was performed using default parameters. Additional information about Cell Ranger can be found at https://support.10xgenomics.com/single-cell-gene-expression/software. Upon completion, Cell Ranger produced an output directory per file that contains the following: bam file (binary alignment file), HDF5 file (Hierarchical Data Format) with all reads, HDF file containing just the filtered reads, summary report (html and csv), cloupe.cloupe (a file for the 10x Loupe visual browser). HDF5 files (filtered) were then uploaded into the R statistical programming language (version 3.6.0) using the Seurat package (version 3.0) where normalization, scaling, integration and reference-based label transfer was performed for cell typing/classification.
Individual HDF5 files were loaded into the R statistical programming language (version 3.6.0) using Bioconductor (version 3.1.0) and the Seurat package (version 3.1.5). For simplicity, sample names were captured as a list in R and iteratively processed within a loop (refer to https://satijalab.org/seurat/ for more information). Within the loop, samples were normalized with the NormalizeData function followed by the FindVariableFeatures function with parameters: vst selection method and 2000 features. Label transfer was performed using previously published procedures (Stuart et al., 2019) . Data was labeled with the Seurat 2 3 reference dataset and checked for expected DEG by cell type (Supplemental Table S7 ). Labeling included the FindTransferAnchors and TransferData functions performed in the Seurat package. Labeling was performed in specified sequence where information acquired from the previous labeling was accumulated then set as anchors for the next sample, looping by time point across patients. The time point to time point variation first observed within the Nanostring data was used as the basis for this labeling strategy.
After label transfer was complete, we calculated read depth, mitochondrial percentage, and number of UMIs per sample. Despite normalization being performed during the label transfer process, the raw counts were stored with cell labels. After labeling and metrics were recorded, samples were merged together in Seurat using the merge function. The merged data structure was normalized (using NormalizeData and FindVariableFeatures functions) and then saved as an RDS for further analysis. For confirmation of cell types, we used the FindMarkers function (Seurat) and the MAST package to identify DEGs for each identified cell type against all other cell types. Unsupervised cell clustering was also performed separately for Experiment 1 and Experiment 2 using the Louvain / KNN-based method in the Seurat package.
Differential expression analysis from later timepoints compared to two hours was conducted on a per reference-based cell type basis and for selected timepoints on a per Louvain-cluster basis using the MAST package (https://pubmed.ncbi.nlm.nih.gov/26653891/) and the FindMarkers function from the Seurat package (version 3.1.5). Because they were sequenced in different batches, differential expression was run separately for Experiment 1 Healthy Donors, Experiment 1 SLE Donors, and Experiment 2 Healthy Donors.
Pathway analysis was performed on single-cell transcriptomics data from Experiment 2. Individual cells were labeled by cell type and clustered as indicated in Methods/Single-Cell Transcriptomics/Cell Classification. As an example of pathway analysis, two clusters exhibiting extreme and inverse time-dependent bias (e.g. present in a cell type at 2 hours but not present for the same cell type at 18 hours, and vice versa) were selected and the corresponding data was submitted to Ingenuity Pathway Analysis (Qiagen, release: Summer 2020) using log2 fold change values and adjusted p-values calculated by differential gene expression analysis. The clusters selected for analysis were (colors as shown in Figure 5B ): magenta ("early") and teal ("late"). Similar analysis was performed for CD4+ memory T cells: pink ("early") and brown ("late").
centrifugation and processing 