Multimodal single-cell analysis of cutaneous T-cell lymphoma reveals distinct subclonal tissue-dependent signatures

Insights into the molecular progression of cutaneous T-cell lymphoma (CTCL) have revealed a significant degree of heterogeneity in affected tissues between different patients and disease stages and prospective intraindividual variability.^1-10 In addition, single-cell profiling of advanced stage tumor biopsies and circulating Sézary syndrome (SS) cells has further highlighted substantial inter- and intratumoral variability.^11,12 While these approaches have provided important insights into the genomic instability and the transcriptional heterogeneity of CTCL, there is a gap in the understanding of disease progression between skin and blood, and little is known about the degree to which malignant T cells in the blood resemble malignant T cells in lesional skin within leukemic disease.^13-15 In the present study, we interrogate multiple modalities of information in single cells from matched skin and blood samples of patients with leukemic disease and healthy controls within a single workflow. By using expanded CRISPR-compatible cellular indexing of transcriptomes and epitopes by sequencing (ECCITE-seq),¹⁶ we sought to compare the molecular profile of malignant clones residing in the skin and circulation of these patients across gene and protein expression modalities. We apply inferred copy number variation (CNV ) and phylogenetic analysis to examine subclonal heterogeneity to gain insights into the evolution and the relationship of malignant clones across tissues.

Patient and control samples were collected at New York University Langone Health Medical Center, Columbia University Irving Medical Center, and Roger Williams Medical Center in accordance with protocols approved by the New York University School of Medicine Institutional Review Board and Bellevue Facility Research Review Committee (IRB#i15-01162), Columbia University Irving Medical Center Institutional Review Board (IRB#AAAQ8751), Roger William Medical Center Protocol Review Committee (IRB#15-520-99), and Rhode Island Hospital Institutional Review Board (IRB#4159-16). Patients with CTCL were diagnosed according to ISCL/EORTC classification criteria, and all patients included in the study had detectable leukemic involvement and staging of mycosis fungoides (MF) stage IV or SS. Most patients included were not undergoing any topical or systemic therapies for at least 2 weeks before the time of biospecimen collection. A summary of the clinical notes, including past therapies for each patient, is included in supplemental Table 1, available on the Blood Web site. After written informed consent was obtained, skin biopsy and peripheral blood samples were harvested. All authors had access to de-identified patient clinical data.

Peripheral blood mononuclear cells (PBMCs) were isolated from the blood of patients and healthy controls by gradient centrifugation using Ficoll-Paque PLUS (GE Healthcare) and Sepmate-50 tubes (Stemcell). Buffy coat PBMCs were collected and washed twice with phosphate-buffered saline (PBS) 2% fetal bovine serum (FBS). Skin biopsies were digested in 0.25 mg/mL of Liberase TL and 0.5 mg/mL DNase (Sigma) diluted in RPMI and incubated for 1.5 hours at 37°C. Skin dissociated cells were filtered through a 70-μM cell strainer. All cells from single-cell suspensions were cryopreserved in freezing medium (40% RPMI 1640, 50% FBS and 10% dimethyl sulfoxide). Cryopreserved PBMCs and skin-dissociated cells were thawed for 1 to 2 minutes in a 37°C water bath, washed twice in warm PBS 2% FBS, and resuspended in complete media (RPMI 1640 supplemented with 10% FBS and 2 mM L-Glut) on the day of the experiments. Biopsies from 2 patients did not yield sufficient T cells to be included in analysis.

Antibodies used for cell hashing and for our ECCITE-seq antibody-derived tags (ADT) panel were covalently and irreversibly conjugated to barcode oligos by iEDDA-click chemistry (“home conjugated”) as previously described^16-18 or purchased as TotalSeq-C reagents (BioLegend). Concentrations and clones used for all conjugated antibodies for staining are listed in supplemental Table 4.

Cells were stained with barcoded antibodies and loaded onto the 10x Chromium single-cell workflow (10X Genomics) as previously described¹⁶ with a few modifications. Approximately 1 million to 3 million cells per sample were resuspended in CITE-seq staining buffer (2% BSA, 0.01% Tween in PBS) and incubated for 10 minutes with Fc receptor block (TruStain FcX, BioLegend and FcR blocking reagent, Miltenyi) to prevent antibody binding to Fc receptors. Cells were then incubated with hashing and ECCITE-seq surface panel antibodies for 30 minutes at 4°C. After staining, cells were washed 3 times in PBS containing 2% BSA and 0.01% Tween, followed by centrifugation (300 ×g for 5 minutes at 4°C) and supernatant aspiration. After the final wash, cells were resuspended in PBS and filtered through 40-µm cell strainers.

Stained cells from each sample were pooled and loaded into the 10x Chromium Single Cell Immune Profiling workflow according to the manufacturer’s instructions with a few modifications as previously described.¹⁶ Additional modifications included 1. Addition of 2 pmol of TotalSeq-C BioLegend protein-tag additive (CTCGTGGGCTCGGAGATGTGTATAAGAGACG) during cDNA amplification PCR for pools stained with TotalSeq-C antibodies and 2. N7NY_xx indexing primers, as described by Mimitou et al,¹⁶ were used along with SI-PCR primer to make libraries for TotalSeq-C antibodies. Libraries were pooled and sequenced on Illumina NovaSeq 6000, MiSeq, or HiSeq 2500.

FASTQ files from the 10x libraries were processed using the count module of Cell Ranger pipeline, version 3.0.1 (10x Genomics). FASTQ files from the T-cell receptor (TCR) α/β and TCR γ/δ libraries were processed using the vdj module of Cell Ranger pipeline, version 3.0.1 (10x Genomics). For ADT libraries and hashtag oligos (HTO) libraries, the tag quantification pipeline version 1.3.2, available at https://github.com/Hoohm/CITE-seq-Count, was used to tally the counts from reads. A Hamming distance of 1 was used as criteria for read sequence matches to preassigned barcode sequences. All reads from the 5' mRNA, TCR α/β, and TCR γ/δ libraries were aligned to the GRCh38 reference. Raw UMI gene expression (GEX), TCR α/β, TCR γ/δ clonotype, ADT, and HTO count matrices were obtained. Samples were demultiplexed by their HTO using Seurat HTODemux with positive quantile parameter set to 0.999.

Scrublet¹⁹ was used to detect and remove intrasample cell duplicates. We parameterized the model using a 7% expected doublet rate, 30 principal components (default), and a minimum gene variability of the 85th percentile (default). We performed manual thresholding by using a doublet score of 0.2.

The count matrices of all modalities were loaded into a Seurat v3 object (Seurat version 3.0.2). For GEX data, we normalized the count data using LogNormalize and scaled the resulting normalized values using ScaleData. For all other modalities (ADT, HTO, TCR α/β, TCR γ/δ), we normalized the count data using centered-log ratio (margin = 2) and scaled the resulting normalized values using ScaleData.

We applied weighted-nearest neighbor (WNN) analysis on our ECCITE-seq data, enabling the integrative analysis of GEX and ADT modalities in the same cells as previously described.²⁰ k-Nearest neighbors graphs (k = 20) were first constructed, and the neighbors of each cell in each modality were identified. For graph construction, the first 30 principal components (PCs) of the log-normalized values from the GEX modality and the first 10 PCs of the centered-log ratio values from the ADT modality were used. The dimensionally reduced molecular profiles of the k-nearest neighbors of each cell were then averaged and compared with the actual values to obtain the residuals. This procedure was done in a within-modality (RNA-RNA, protein-protein) and cross-modality (RNA-protein, protein-RNA) manner. Next, for each of these residuals, they were converted into affinity-based similarities using Uniform Manifold Approximation and Projection (UMAPs) locally adaptive kernel. Finally, the ratios of the within-modality and cross-modality affinities were softmax transformed to obtain the cell-specific modality weights. The weighted similarity then formed the basis of the construction of the WNN graph. Using the integrated WNN graph, we constructed t-SNE visualizations and performed modularity clustering from it.

TCR α/β reads attributed to each cell barcode were grouped, assembled into a single contig, and then had their V, D, and J segments annotated. Only clonotypes with productive and full-length contigs were retained for analysis. TCR α/β and TCR γ/δ clonotype count matrices were obtained and appended alongside GEX and ADT modalities for simultaneous analysis.

InferCNV version 1.2.1 was applied for analysis of individual patient samples using malignant cells (defined by their expanded TCR β clonotype, distinctive clustering, and GEX patterns) (supplemental Figures 3 and 4) from both the blood and skin as the targets, and nonmalignant CD4+ T cells were used as the reference. For subcluster analysis, the analysis_mode parameter was set to “subclusters.” Default parameters (cutoff = 0.1, cluster_by_groups = TRUE, denoise = TRUE, HMM = TRUE) were used for these analyses.

Using the de-noised CNV inference hidden Markov model results, we built copy number profiles for each patient at 500-Kb resolution. Neighbor-joining trees of malignant T cells were then constructed and rooted with an individual nonmalignant T cell using R package ape version 5.3. Bootstrap analysis of the NJ tree was performed using 100 replicates to estimate the confidence of the clades. Visualization of the tree as a phylogram was done using R package ggtree version 2.0.1.

Pseudotime trajectory analysis was performed on each patient with CTCL with matched blood and skin samples (SS1 to SS4 and MF stage IV) using Monocle 2 version 2.14.²¹ Briefly, moderately expressed genes identified using detectGenes (min_expr = 0.1) and at least 20 cells expressing the genes were used for the analysis. PCA was performed, and the top 20 components were used for t-SNE visualization. Ordering genes were defined based on DifferentialGeneTest on the tissue of origin and used for trajectory construction (DDRTree)^21,22 and pseudotime analysis. We rooted the state based on the CNV subtype that has the shortest branch length as compared with nonmalignant CD4+ T cells.

Two gene sets—resting CD4+ T cells and activated CD4+ T cells upon CD3/CD28 stimulation in human tissues—were retrieved from Szabo et al.²³ We used Seurat function ScoreModules to score each malignant cell from both blood and skin on each gene set.

Code for the pipeline to pre-process ECCITE-seq data are available on GitHub: http://github.com/ouyang-lab/CTCL. All downstream analysis will be made available on request.

We used ECCITE-seq to examine clonal malignant T cells and survey other leukocytes in the blood and skin of 5 patients with leukemic MF and SS and in the blood of 2 additional patients with SS (supplemental Table 1). While sequencing-based capture of surface epitopes recapitulates distribution of markers used for fluorescence -activated cell sorting analysis of malignant T cells in CTCL (supplemental Figure 1), ECCITE-seq enables a far more detailed analysis of the samples with information on >50 surface epitopes, single cell transcriptome and V(D)J repertoire.¹⁶ As expected, we detected clonally expanded malignant T cells only in CTCL patient-derived specimens (Figure 1A-B; supplemental Figures 2-4). Differences in malignant burden and in the distribution of other leukocyte populations were observed between individual patients (Figure 1B; supplemental Figure 2; supplemental Tables 2 and 3). The malignant T cells within each patient were defined based on their expanded TCR β CDR3 and on distinct GEX and ADT profiles clustering them away from polyclonal nonmalignant T cells in both tissues (Figure 1A; supplemental Figure 2). Malignant GEX profiles included dysregulated expression of hallmark genes frequently affected in leukemic MF and SS, including SATB1 , DPP4 (CD26), CD7, TOX, PSL3, HDAC9, DNM3, CDO1, KIR3DL2, and GATA3 (supplemental Figures 3 and 4). Additionally, the malignant T cells in all patients with CTCL with matched lesional skin and blood shared identical TCR clonotypes between tissues, indicating distinct clonal involvement across tissues in late-stage CTCL.

We then reclustered our CTCL samples, only including the malignant T cells, and focused on their transcriptional profiles. Our analysis detected transcriptional heterogeneity at the single-cell level among malignant T cells from the same patient, allowing for grouping into distinct clusters (Figure 1C-D; supplemental Figure 5). Interestingly, when malignant T cells from skin and blood were analyzed together, most of the subclonal transcriptional heterogeneity was observed in the blood, whereas the skin-derived malignant T cells mainly constituted a single subcluster (Figure 1C-D; supplemental Figure 5). This was consistent across all patients investigated; together, these results exposed a high degree of transcriptional heterogeneity in blood-derived malignant T cells of multiple patients with CTCL and a unique transcriptional bias generally grouping skin-derived malignant T cells in their own subcluster.

SS is known to encompass a significant degree of genetic instability across patients, portrayed by an abundance of CNVs.^1-32425 To determine if the intraclonal malignant transcriptional heterogeneity in our CTCL cohort reflects the presence of distinct malignant subclones, we inferred CNVs at the single-cell level in clonal malignant T cells from single -cell RNA-sequencing (scRNA-seq) data as has been previously demonstrated in other tumor types²⁶ (Figure 2A-C; supplemental Figure 5). This technique (InferCNV) was previously used in the analysis of uveal melanoma to ascertain clonal structure in tumors and reveal noncanonical large-scale CNVs with potential implications in tumor progression.²⁷ CNV inference methods use GEX data to find broad, contiguous sections of consistent copy number changes. We took advantage of this method to identify broad CNV events to distinguish major subclonal malignant populations. We corroborated our results using an alternative method of CNV profiling from scRNA-seq (copy number karyotyping of aneuploid tumors)²⁸ (supplemental Figure 6). One limitation of these algorithms involves possible false discovery of CNV patterns if significant GEX differences occur between the chosen reference and target cells in the absence of genetic aberrations. We reinforced our confidence in inferring CNV events from scRNA-seq data by applying the method on scRNA-seq from a patient with SS not otherwise included in the study and comparing with matched whole-exome sequencing and methylation microarray datasets. This comparison revealed high concordance in the CNVs detected by scRNA-seq, whole-exome sequencing, and methylation microarray, albeit with low sensitivity for small-scale aberrations, as would be expected (supplemental Figure 7). To independently validate our inferred CNVs, we also measured loss of heterozygosity stemming from CNV events in our individual patients with CTCL (supplemental Figure 8). Our analysis of skin and blood revealed multiple subclonal populations within each patient’s malignant T cells, defined by varying degrees of acquired CNVs in both tissues (Figure 2A-C). Each patient also contained a unique combination of CNVs within their malignant subclones, reinforcing the wide range of genetic aberrations found in late-stage CTCL tissues. Similar CNV profiles are often shared between individual malignant subclones across the tissues, suggesting that genetic heterogeneity alone does not dictate tissue localization of malignant cells. Notably, the CNV defined subclones in 5 out of 7 patients (SS1-SS4 and SS6) showed a moderate degree of concordance (adjusted Rand index: 0.28-0.50) with the previously transcriptionally defined clonal subclusters (Figure 2D; supplemental Figure 5) supporting InferCNV’s predictions. Although the other 2 patients (MF stage IV1 and SS6) showed a lower degree of alignment (adjusted Rand index: 0.08-0.10), most cells from each defined state still visually aligned with each other.

Because CTCL is thought to arise from skin effector/memory T cells and progress to leukemic involvement in a subgroup of patients,^13,29,30 we sought to examine whether the evolutionary relationship between CNV-defined subclones in each patient can identify the tissue of origin of malignant cells and provide a map of clonal evolution in leukemic disease. To this end, we applied a phylogenetic approach in which we constructed neighbor-joining trees as well as minimum evolution and average linkage distance–based phylograms for each patient’s malignant T cells, rooted by the polyclonal nonmalignant T cells from the same patient (Figure 2C; supplemental Figures 5, 9, and 10). We observed that in most of the matched samples, there was no clear monodirectional phylogenetic relationship between skin- and blood-derived subclones to suggest a tissue of origin for disease, with distinct trajectories of clonal evolution among individual patients. Our analysis showed branches of clonal evolution leading to advanced skin and blood subclones in each patient. Advanced subclones also shared several common intermediate-late branching points regardless of tissue distribution, indicating that closely related subclones could reside in either affected tissue. Interestingly, only the patient diagnosed with MF stage IV1, who did not meet clinical criteria for SS diagnosis, showed a more clear distinction in subclonal trajectories between tissues. This complex and variable phylogenetic landscape based on CNV accumulations reveals a snapshot of the development of multiple related but unique subclones in the skin and blood of patients with SS, with both tissues frequently sharing similar branches of evolution. Collectively, these data highlight the highly heterogeneous nature of late-stage CTCL at the subclonal level and highlight a clinically relevant challenge presented by the emergence of distinct subclonal populations across multiple tissues.

In contrast to their blood-derived counterparts, skin-derived malignant T cells from our patients with CTCL showed a strong tendency to group either as a single main subcluster (SS1, SS3-MF stage IV1) or as a separate island of cells away from circulating cells (SS2) when analyzed by their transcriptional profile (Figure 1C; supplemental Figure 5). To assess the transcriptional relationship between malignant skin- and blood-derived T cells from these patients, we applied a transcriptional trajectory analysis of matched samples as previously described²¹ (Figure 3A,C; supplemental Figure 11). This analysis revealed transcriptional trajectories that generally started with subclones having the fewest CNVs and ended with subclones containing more extensive CNVs. We also observed a strong transcriptional bias, with skewing of all skin-derived subclones, toward the end of the branches that started at blood-derived subclone trajectories. This is in contrast with our previous observations using phylogenetic analysis of the inferred CNV subclones, in which we found an intermixed distribution among skin and blood subclonal branches in regards to their path of evolution and their common ancestors (Figure 2B; supplemental Figure 5). These results suggested the presence of a strong tissue-dependent transcriptional bias in malignant T cells from the skin, effectively masking subclonal differences.

To determine the underlying mechanism behind this transcriptional bias observed in skin-derived malignant T cells, we performed paired transcriptional comparisons between matched skin- and blood-derived malignant T cells from 5 of our patients with CTCL (Figure 3B-C; supplemental Figure 12). This comparison revealed a consistent and strong upregulation of several T cell activation-, TCR ligation-, and mitogen-induced transcripts as well as regulators of cell cycle in malignant T cells from skin as compared with their blood-derived counterparts. These genes included PDCD1 (PD-1), NR4A1, NR4A2, NR4A3, TNFRSF18 (GITR), TNFRSF4 (OX-40R), BRD2, REL, DUSP2, DUSP4, GEM, and CDKN1A among others. Additionally, upregulated transcripts in blood-derived malignant T cells compared with their skin-derived counterparts included KLF2, an important factor in maintaining T-cell quiescence,³¹ and TCF7 and SELL (CD62L), known markers of resting T cells.^32,33 Because PD-1 and CD62L were included in our ADT surface protein panel (supplemental Table 4), we examined whether the tissue-dependent expression of these proteins was also evident in individual malignant T cells (supplemental Figure 13). Our results revealed significantly higher levels of surface PD-1 expression in skin-derived malignant T cells compared with blood counterparts in 4 of 5 patients with CTCL (SS1-4) and significantly higher levels of CD62L expression in blood-derived malignant T cells compared with their skin counterparts in 4 of 5 patients with leukemic CTCL (SS2-5) in line with tissue-specific distribution of activation signature.

The abundance of differentially expressed transcripts relating to differences in activation status between malignant T cells depending on their tissue location prompted us to investigate whether independent “T-cell activation status” gene sets were enriched in malignant T cells from skin and blood (Figure 3A,D; supplemental Figure 11).²³ Notably, skin-derived malignant T cells had a consistently higher T-cell activation score compared with blood-derived malignant T cells, whereas the opposite was true for the resting T-cell score. It resulted in a markedly increased activation-to-resting-score ratio in malignant T cells from skin as compared with their blood-derived counterparts in each patient (Figure 3D). Importantly, this activated signature is unique to the lesional skin microenvironment of CTCL, because there was no difference in the activation to resting score ratio between skin-derived CD4+ T cells from the nonlesional skin of a healthy control and a patient with psoriasis compared with their blood counterparts (supplemental Figure 14).

To explore whether this chronically activated phenotype of malignant T cells in the skin tumor microenvironment coincides with a high proliferative potential, we applied a cell cycle scoring analysis using canonical markers of G1, S, and G2/M as previously described³⁴ and superimposed them on the transcriptional trajectories of malignant T cells (supplemental Figure 15). Our results revealed that the frequency of highly proliferative malignant T cells in the skin of patients with CTCL was higher than that found in circulation. These results are consistent with observed activation of malignant T cells in the skin and indicate that the skin tumor microenvironment provides signals conducive to frequent malignant T-cell activation and expansion. Of note, regressing out tissue-dependent transcriptional signatures allowed malignant subclones, based on inferred CNV profiles from skin, to cluster together with their most phylogenetically closely related subclones in blood (supplemental Figure 16). Overall, these results show consistent transcriptional signatures delineating malignant subclones from skin and blood, highlighting the importance of the tumor microenvironment in shaping the transcriptional programs that contribute to distinct functional responses.

SS is known to be a heterogeneous disease, with different patients demonstrating highly variable transcriptional profiles,^6,11,35-37 but not much is known about the relationship between malignant T cells from the skin and blood in these patients. Our findings reveal the presence of a distinct T-cell clone identified in both the matched skin and blood of patients with MF stage IV and patients with SS with multiple CNV-defined subclonal populations in both tissues. In accordance with previous studies,^1,3,38 many of these subclones displayed recurrent CNVs, particularly involving gains in 8q (containing MYC , 5 of 7 patients) and 17q (5 of 7 patients) and losses in 10q (containing PTEN, 4 of 7 patients) and 17p (containing TP53, 5 of 7 patients) that could be found in both tissues. This broad assortment of CNVs was accompanied by a high degree of transcriptional heterogeneity between clonal malignant T cells across patients (supplemental Figure 17). The presence of a highly branching subclonal evolution regardless of tissue was also reflected in our phylogenetic analysis, in which we saw no clear linear relationship between skin- and blood-derived subclones. GEX analysis using T-helper (Th) subtype differentiation markers between tissues also revealed preferential expression of Th2 and Th17 signature transcripts across malignant T cells, consistent with previous reports on Th identity of malignant cells in CTCL^39-41 (supplemental Figure 18). These results highlight a close genetic relationship between subclonal populations in circulation and skin, frequently sharing branches of clonal evolution, and suggest that migration between tissues is continuous, not a one-time event. While one of our patients (MF stage IV1) showed clear phylogenetically related subclones based on tissue distribution, that individual was the only patient not considered to have SS based on lack of clinical criteria, suggesting that recirculation of evolved malignant clones may be a characteristic of SS. This highly diverse landscape at the subclonal level provides insight regarding the refractory nature of advanced-stage CTCL,⁴² in which malignant subpopulations can be resistant to treatment,⁶ and highlights the need for a personalized, multipronged therapeutic approach wherein the complete landscape of subclonal diversity in each patient should be considered. Because most patients with SS are treated with sequential single-agent therapies,³⁰ our data provide additional rationale behind a polytherapeutic approach.

Our differential GEX analysis revealed robust differences between malignant T cells from blood and skin and highlighted marked enrichment of transcripts related to activation status and cell cycle progression. This indicates that the tumor microenvironment in CTCL can dictate the transcriptional program of malignant T cells. Several lines of evidence have pointed to the critical role of the skin microenvironment in the pathogenesis of CTCL, a site of interaction between skin pathobionts and immune cells and as a site of chronic activation and proinflammatory signaling that promotes the selection of malignant CD4+ T cells.^11,43-50 A report by Cristofoletti et al⁵¹ also proposes that skin-derived T cells in SS have a higher proliferative index compared with their blood-derived counterparts. Together with the highly consistent differences we observe in expression of activation- and mitogen-induced transcripts in malignant T cells found in the skin and blood of patients, these findings indicate that the skin microenvironment of CTCL promotes a transcriptional landscape driving sustained activation that can lead to rapid malignant expansion.

Because treatment failure is a common occurrence in patients with SS,⁴² it is imperative to improve our understanding of the mechanisms that promote disease heterogeneity, specifically the emergence of multiple malignant subclones. Our results suggest that the skin microenvironment can greatly influence the transcriptional state of malignant T cells and may promote evolution of malignant subclones that can readily disseminate through the circulation. Our study provides a foundation for further understanding the molecular drivers that orchestrate tissue distribution in the skin and blood of patients with advanced-stage CTCL at the single-cell level and makes the case for extending this multimodal approach to interrogate the key drivers of both early skin-confined and progressive disease.

The authors thank New York University Genome Technology Center for technical assistance and support.

S.B.K. was supported by the National Institutes of Health, National Heart, Lung, and Blood Institute R01 grant (HL-125816), LEO Foundation grant (LF-OC-20-000351), New York University Cancer Center pilot grant (P30CA016087), the Judith and Stewart Colton Center for Autoimmunity pilot grant, and a grant from the Drs. Martin and Dorothy Spatz Foundation. A.S. was supported by HHMI Medical Scholar Fellowship. P.S. was supported by the National Institutes of Health, National Human Genome Research Institute (R21HG009748) and the Chan Zuckerberg Initiative (HCA-A-1704-01895). T.B.B. and N.Ø. were supported by the Danish Cancer Society (Kræftens Bekæmpelse), the Danish Council for Independent Research (Danmarks Frie Forskningsfond), and the LEO Foundation. M.E.K. was supported by a grant from the Drs. Martin and Dorothy Spatz Foundation.

Contribution: A. Herrera and S.B.K. initiated the project and were responsible for experiment design with input from E.P.M., N.Ø., A. Heguy, and P.S., E.P.M. and A. Herrera performed the majority of experiments; A.S., J.U.S., L.J.G., M.E.K., D.G., M.B., K.H., and J.L. were directly involved in patient recruitment, sample collection and study design; A.C. performed the bioinformatics work under the supervision of Z.O. with input from A. Herrera, K.V.R., T.B.B., and S.B.K; and A. Herrera, A.C., T.B.B., and S.B.K. wrote the paper with all authors providing input on the manuscript.

Conflict-of-interest disclosure: P.S. is listed a co-inventor on a patent application related to ECCITE-seq (US provisional patent application 62/515-180). The remaining authors declare no competing financial interests.

Correspondence: Sergei B. Koralov, Department of Pathology, New York University School of Medicine, 550 First Ave, MSB 531A, New York, NY 10016; e-mail: sergei.koralov@nyumc.org; and Terkild B. Buus, Department of Immunology and Microbiology, University of Copenhagen, Blegdamsvej 3B, Bldg. 07-12-62, DK-2200 Copenhagen, Denmark; e-mail: terkild.buus@sund.ku.dk.