Clonal Expansion and Interrelatedness of Distinct B-Lineage Compartments in Multiple Myeloma Bone Marrow

Multiple myeloma is a malignancy of the B lineage characterized by the accumulation of clonal plasma cells in the bone marrow. Most, if not all, cases develop from a precancer called monoclonal gammopathy of undetermined significance (MGUS; refs. 1, 2). Due to the presence of identical somatic mutations of the immunoglobulin genes and completed class-switch recombination, the disease is assumed to originate at the germinal center or post-germinal center stage of B-cell development (3–5). Current thinking holds that the origin of this malignancy includes alterations in the bone marrow microenvironment (6, 7), vasculogenesis (8), the T-, B-, NK-, and dendritic cell compartments (9–11), as well as genetic aberrations within the malignant cells themselves (12–16), resulting in variable degrees of intraclonal heterogeneity.

Intraclonal heterogeneity usually refers to accumulating genetic and epigenetic aberrations that drive the progression from MGUS to symptomatic multiple myeloma (17, 18). Clonal heterogeneity has also been detected at the immune phenotype level and certain phenotypes have been shown to be partly associated with genetic aberrations, chemotherapy-resistance, and proliferation advantages (16). Nevertheless, identification of phenotypic diversity has been challenging. Although flow cytometry can detect high-dimensional immune phenotypes on the single-cell level, clonality has usually been determined by immunofixation or immunoglobulin sequencing in bulk populations resulting in the loss of single-cell resolution and accuracy.

Neither RNA- nor DNA-based bulk immunoglobulin gene sequencing can quantify clone frequencies or assigned immune phenotypes to the identified clones could be inaccurate as a result of contamination with other cells (19). Yet the success of B cell-targeting cellular therapeutics in subsets of patients (20, 21) suggests the involvement of cells phenotypically different from the majority of aberrant-phenotype plasma cells in multiple myeloma pathophysiology.

Knowledge of clonal diversity and phenotype is critical to understanding the pathogenesis of multiple myeloma and may aid identification of targets for therapeutic intervention. We combined phenotypic tracking of B-lineage clones using single-cell multiparameter FACS index sorting with single-cell next-generation sequencing (NGS) of immunoglobulin genes in order to assess the phenotypic diversity of any B-lineage clone of choice in multiple myeloma. Index sorting allows the isolation of hundreds to thousands of single cells without upfront selection of specific immune phenotypes, yet fluorescence intensities for all markers included in the panel are recorded and can be read out for each individual cell retrospectively. We developed a methodology that combines DNA-barcoding and high-throughput NGS of immunoglobulin light chain genes from hundreds of single B-lineage cells with 13-parameter FACS index sorting to phenotypically track individual clones on the single-cell level in human bone marrow.

We find that the immune phenotypes of these molecularly defined, predominant multiple myeloma clones range from aberrant-phenotype multiple myeloma cells to normal-phenotype B-lineage cells. The parallel expansion of a variety of less predominant B-lineage clones follows the same pattern.

The majority of the less predominant B-lineage clones phenotypically clustered with normal-phenotype B cells and showed low somatic mutation rates in the light chain genes and silent nucleotide differences within the CDR3, suggesting that their expansion is driven by antigen-dependent physiologic rather than malignant processes.

Fresh heparin-anticoagulated bone marrow from three multiple myeloma patients (MM1-3; Table 1; Supplementary Table S1) and peripheral blood from one healthy volunteer were acquired after they gave written informed consent in accordance with federal and local human subjects regulations (Institutional Review Board ID 25310 to LH). Mononuclear cells were isolated using Ficoll–Paque PLUS (GE Healthcare) centrifugation, resuspended in cell culture medium (Supplementary Table S2) containing 10% dimethyl sulfoxide and 50% fetal bovine serum (both Sigma Aldrich), and cryopreserved.

Cells were thawed and stained for FACS analysis/sorting (Supplementary Table S3) according to the manufacturer's instructions. For the detailed gating strategy, see Supplementary Fig. S1. Plates were preloaded with 1.4× first-strand buffer for Superscript III RT (Invitrogen) and chilled during single-cell index sorting on a FACS Aria II or Aria Fusion cell sorter (BD). Plates were immediately transferred to dry ice after sorting.

Superscript III (Invitrogen) reverse transcription (RT) was done on index-sorted single cells in a final volume of 14 μL/reaction. Octylphenoxypolyethoxyethanol (IGEPAL CA-630, Sigma) was added at 1% final concentration and RT primers (Supplementary Table S4) were added at 257 nmol/L final concentration. Plates were incubated at 65°C for 3 minutes, 25°C for 3 minutes, and placed on ice to add 7.9 mmol/L DTT (Invitrogen), 429 μmol/L dNTPs (Thermo Scientific; all final concentrations), and 0.3 μL Superscript III RT. RT was done at 37°C for 1 hour and 70°C for 15 minutes. First and second amplification reactions were done using Phusion High-Fidelity DNA Polymerase (NEB). Final PCR volumes were 15 μL using 4 μL cDNA or 1 μL first amplification products as templates. dNTPs were used at 433 μmol/L and DMSO was added at 5% final concentration to the second amplification PCR. First amplification and second amplification (barcoding) primers were added to a final concentration of 333 nmol/L each. Cycling conditions: first amplification: 98°C, 60 seconds, 3× (98°C, 45 seconds; 45°C, 45 seconds; 72°C, 105 seconds), 30× (98°C, 30 seconds; 50°C, 45 seconds; 72°C, 90 seconds), 72°C, 7 minutes; second amplification: 98°C, 30 seconds; 30× (98°C, 30 seconds; 50°C, 45 seconds; 72°C, 45 seconds), 72°C, 5 minutes. Primers for RT, first, and second amplification were adapted from Wang and colleagues (22) and modified for molecular barcoding and Illumina sequencing. The barcoding primers add a nucleotide sequence specific for plate, column, and row to each single amplified transcript (Supplementary Table S4). Barcodes were designed with differences in at least two nucleotides when compared with all other barcodes in the panel.

Third amplification reactions were done in 15 μL final volume using a multiplex PCR kit (Qiagen) and 1 μL of second amplification product as template. Paired-end primers adding the remaining sequences of the Illumina adapters were used at 400 nmol/L final concentration and cycling conditions were: 95°C, 15 minutes; 28×(95°C, 30 seconds; 60°C, 45 seconds; 72°C, 90 seconds), 72°C, 10 minutes.

After the third amplification, PCR products from all wells were combined, gel purified, and sequenced on an Illumina MiSeq (Illumina) using Illumina V3 kits (Illumina).

Illumina sequencing reads were assembled using PEAR (23) and a cut off of at least 200 nucleotides for assembled reads. Python scripts identified barcodes in every single read (only perfect matches were accepted), filtered sequences by aligning against germline-encoded V genes using bowtie2 (24), built consensus sequences with MUSCLE (25) to upload them to IMGT/HighV-QUEST (26) for V, J gene assignment, detection of somatic mutations, and CDR3 identification. Read number cutoffs were determined based on background reads in empty wells in analogy to our methodology for T-cell receptor sequencing (ref. 27; Supplementary Fig. S2). Sequence alignments were done running MUSCLE out of MEGA (28). Nucleotide sequences of individual cells of the predominant clones of all three patients have been made publicly available in GenBank, accession number: KY177746-KY178269.

Sequence logos were created with WebLogo 3 (29).

For t-SNE (30) analyses, compensated FACS data and index sort information were exported from FlowJo (TreeStar) and DiVA (BD) software and combined using R (31). The Barnes-Hut implementation of t-SNE was run using the arguments “–v –d 2 –p 10” on combined FACS data exported from gates shown in Supplementary Fig. S1D and index sort data. t-SNE visualization was done using the ggplot2 package in R. For the visualization of individual marker expressions, fluorescence values were trimmed at the 1st and 99th percentiles, scaled, log₁₀-transformed, and color coded with yellow representing the lowest and red representing the highest expression levels.

For two-dimensional FACS plots of single clonal cells, FACS data from the whole sort dataset (contours in gray) were overlaid with fluorescence data of single clonal cells (filled black circles) of the indicated clones in R. In two-dimensional FACS plots of individual clones, negative fluorescence values were set to 0 for visualization purposes.

In multiple myeloma, flow cytometry can identify malignant cells through their aberrant antigen expression profiles (16, 32, 33). Although a variety of markers characterize this disease, variation is considerable between and within individual patients. To better capture this variation, we used a nonlinear principal component analysis algorithm, t-SNE (30) for the visualization of multidimensional (13 markers and forward-/side-scatter characteristics; Supplementary Table S3) flow cytometry data of multiple myeloma bone marrow cells from three heterogeneous patients (MM1-3; Fig. 1 and Table 1). Single data points in the t-SNE maps represent single cells. Coordinates were calculated taking into account all pairwise distances in the multidimensional datasets. Malignant plasma cells in MM1 were CD38^hiCD27⁺ with aberrant coexpression of CD56, but lacked CD45 and CD20 expression. Plasma cells in MM2 were CD38^hi but did not coexpress CD27, CD56, CD45, or CD20. The aberrantly differentiated cells in MM3 were CD38^hiCD27^loCD56^lo. In each of the three multiple myeloma cases the plasma cell populations occupied multiple compartments in the t-SNE maps, underlining their phenotypic heterogeneity even within individual patients (Fig. 1). MM3 received treatment until two months before the bone marrow aspiration done for this study.

Current flow cytometry approaches to multiple myeloma can identify diversity in the malignant cell populations through their aberrant marker profiles, even in the absence of information about the individual cells' immunoglobulin sequences. However, these approaches do not allow visualization of the phenotypic distribution of single cells belonging to a particular, molecularly defined B-lineage clone.

Phenotypes and clonal relatedness of individual cells can be identified at the single-cell level. Therefore, we combined single-cell FACS index sorting with DNA-barcoding and high-throughput single-cell immunoglobulin light chain sequencing to determine the clonal relatedness of individual B-lineage cells and to define their position on high-dimensional t-SNE phenotypic maps (Fig. 2). We applied a 13-marker FACS panel to bone marrow cells from three multiple myeloma patients and index sorted single B-lineage cells identified by CD19 and CD38 expression characteristics (Fig. 2A, left) into multiwell plates. Index sorting allowed assignment of fluorescence intensities for each fluorchrome in the FACS panel and forward/side scatter characteristics to every cell sorted independently from the sorting gates. In addition to other B-lineage cells, we sorted phenotypically defined multiple myeloma plasma cells (Supplementary Fig. S1) to identify the multiple myeloma immunoglobulin rearrangement. Light chain RNA of each single sorted cell was reverse transcribed, amplified, barcoded with plate- and well-specific DNA barcodes, and sequenced using the Illumina MiSeq system. For the PCR amplifications, we used a nested PCR approach (Fig. 2B) resulting in the amplification of κ and λ light chain genes spanning from the framework region 1 to the constant region. Cells showing identical VJ junction amino acid sequences were considered clonal. Clonality was confirmed by comparing amplified nucleotide sequences (framework region 1 to constant region, approximately 344 nucleotides), including somatic mutations. In parallel, t-SNE maps were generated from the corresponding FACS data and index sort technology allowed the precise mapping of each cell (Fig. 2A, right, clonal cells of the predominant clone were colored in black). Our sequencing approach led to the identification of light chain sequences with a median efficiency of 71% of the sorted cells (range, 32%–84%). Efficiency varied with the stringency of the applied sorting gates. We sought to include all possible B-lineage cells in the sorting gate at the expense of overall sequencing efficiency by allowing selection of a variable percentage of non-B cells (Fig. 2A, left, and Supplementary Fig. S1).

To obtain an experimental estimate of the accuracy of our index sorting, we randomly sorted single T and B cells into 96-well plates and used our sequencing technology to identify B cells that were wrongly assigned a T-cell phenotype. A wrong assignment did not occur in even a single instance and therefore the probability for an index sort mistake was below the experimental detection limit of 0.0045 (Supplementary Fig. S3).

Hundreds of single B-lineage cells from the same multiple myeloma patients whose phenotypes were presented in Fig. 1 were index sorted and sequenced as outlined in Fig. 2. Cells of the predominant clones were overlaid in black filled circles over the corresponding t-SNE plots (Fig. 3A). The monoclonal cells mapped to plasma cell compartments but with aberrant marker expression patterns. In MM1, most cells of the predominant B-lineage clone colocalized with CD38^hiCD27⁺CD56⁺ cells. In MM2, the cells of the predominant clone were CD38^hiCD45⁻ and in MM3 cells of the predominant clone mapped to the CD38^hiCD27^loCD56^lo compartments (Fig. 3A).

In addition to aberrant-phenotype plasma cells, we detected cells that had the same characteristic immunoglobulin nucleotide sequence (frame work region 1 to constant region, including somatic mutations; Supplementary Fig. S4) but mapped to other phenotype compartments. For example, single cells from MM1, when plotted in two dimensions with a selected set of markers (Fig. 3B), showed expression of CD19, CD20, and CD45 alone or in combination without coexpression of CD56. A few single cells of the predominant clones showed phenotypic characteristics similar to normal B rather than multiple myeloma cells in up to 1.2% of all sequenced B-lineage cells (up to 6 cells per patient; for detailed FACS plots of the predominant clones of all multiple myeloma patients see Supplementary Fig. S5). Each single B cell of normal phenotype that belonged to the multiple myeloma clonotype yielded on average 1,137 sequencing reads (range, 89–3228 reads), which was >8-fold more than background (Supplementary Fig. S2). Although multiple myeloma cells have undergone class switch recombination and lack CD20 as well as surface immunoglobulin expression, one B-lineage cell in MM2 shared the characteristic light chain rearrangement and was CD19⁺CD20⁺CD27⁻ in combination with weak surface IgD expression, presenting an immune phenotype similar to naïve B cells (Supplementary Fig. S5).

To confirm the identity of the light chain sequences representing the multiple myeloma clones (predominant clones in Table 2), we sequenced the amplified fragments (framework region 1 to constant region, approximately 344 nucleotides). The sequences, including those nucleotides deriving from somatic mutation (Supplementary Fig. S4), were identical. This suggests that these multiple myeloma clones did not arise in or reenter germinal centers, where variation due to somatic mutation might be expected.

Less predominant clones accompanying the predominant multiple myeloma clones have been identified by immunofixation and on the phenotypic (FACS) or molecular (bulk sequencing) levels (16, 34–36). However, sequencing analyses in bulk are prone to PCR bias and the phenotypic relationships of less predominant clones to the multiple myeloma clones cannot be determined by bulk sequencing. Our methodology allows for quantification of less predominant clones in combination with their immune phenotypes on the single-cell level.

In MM1, we detected a total of 45 expanded B-lineage clones (their CDR3 sequences occurred in two or more sorted cells), 38 in MM2, and 21 in MM3 (Supplementary Table S5). CDR3 sequences and frequencies of the 5 most predominant clones are summarized in Table 2. The majority of the less predominant clones (clones 2–5 in MM1 and MM2, clones 4–5 in MM3) did not map to multiple myeloma plasma cell compartments and showed phenotypic characteristics of normal B cells (CD45⁺CD19⁺CD20⁺CD27^+/−CD56⁻CD138⁻, surface IgD^+/−) as shown for the second predominant clone in MM2 as an example in Fig. 4A. The complete phenotypes of less predominant clones are shown in Supplementary Fig. S6. However, the second and third predominant clones in MM3 showed mixed characteristics of B cells (CD45⁺, partially CD19⁺CD20⁺) and plasma cells (CD38^hi, surface IgD⁻IgG⁻, Fig. 4B, Supplementary Fig. S6). These clones mapped to the same phenotypic compartments that were occupied with the corresponding multiple myeloma cells (see Fig. 3A) and the partial coexpression of CD56 characterized these cells as aberrantly differentiated B-lineage cells (Fig. 4B, detailed phenotypes in Supplementary Fig. S6).

The second most predominant clones in MM2 and MM3 showed identical CDR3 sequences (Table 2). Due to less diversity in immunoglobulin light chain sequences when compared with heavy chains, overlapping CDR3 sequences between individuals is not surprising. In fact, this particular nucleotide sequence of MM2 and MM3 has been previously identified and deposited in the NCBI GenBank database (accession number EU789132). The expansion of this particular light chain in the context of multiple myeloma has to be confirmed in larger cohorts and the underlying cues will be identified in future studies.

The finding of less predominant B-lineage clones in multiple myeloma bone marrow raises the question of whether these clones are the result of a parallel expansion of malignant B-lineage cells or part of a normal antigen-directed response. The light chain gene sequences of all cells of the most predominant clones showed nucleotide sequence identity across the entire VJ gene segment including somatic mutations (see Fig. 5A for numbers of somatic mutations and Supplementary Fig. S4 for sequence logos), indicating their malignant clonal relationship. In contrast, 8 out of the 12 less predominant clones listed in Table 2 showed unique CDR3 sequences and a nucleotide diversity within the CDR3 of individual clonal cells that did not affect amino acid sequences (example in Fig. 5B). Similar patterns could also be observed in peripheral blood B cells from healthy individuals (Supplementary Fig. S7), suggesting that less predominant clones in multiple myeloma bone marrow could arise as part of a normal, antigen-driven expansion (Supplementary Fig. S4 for sequence logos of the five most predominant B-lineage clones from all multiple myeloma samples). Less predominant clones harbored varying (and in 10 out of 12 cases significantly lower, P < 0.05, Fig. 5A) amounts of somatic mutations in their light chain sequences when compared with the predominant clones (Fig. 5A). In line with a lower number of somatic mutations, 7 out of 12 less predominant clones showed surface IgD expression and were CD45⁺CD20⁺, underlining their phenotypic and molecular difference from the predominant clones (shown for clones 1–3 of MM2 as an example in Fig. 5C).

Taken together, the expansion of the most predominant multiple myeloma clones, despite their phenotypic diversity, is part of the malignant monoclonal expansion and shows its phenotypic range. The minor clones in most cases do not show plasma cell phenotypes and seem characteristic of a normal, antigen-driven process.

Estimation of B-cell clonal frequencies and identification of clonal phenotypes in multiple myeloma require the efficient and reliable combination of single-cell technologies.

The application of single-cell methods is especially useful here, as it overcomes the bulk sequencing bias due to the variable number of immunoglobulin gene transcripts per cell. Especially when analyzing bone marrow cells of the entire B lineage, where plasma cells can contain 10–300 times more immunoglobulin RNA than mature B cells (37), bulk sequencing approaches using immunoglobulin mRNA as a template are apt to be especially biased. DNA-based approaches are less affected by varying template copy numbers per cell but are still subject to PCR amplification bias and in general achieve lower efficiencies.

As sequencing efficiency is important for our methodology, we focused on immunoglobulin light chain sequencing, which yields higher efficiencies when compared with heavy chain sequencing. Despite less junctional diversity in light chain than in heavy chain immunoglobulin genes, the substantial amount of somatic mutations in multiple myeloma cells (at average 24 somatic mutations in the most predominant clones in our dataset) allow us to detect clonality (38). Unproductive heavy chain rearrangements can occur in approximately 15% of multiple myeloma patients (39, 40).

The combination of sequencing technology (a median efficiency of 71%) with multicolor (13 parameters) single-cell FACS index-sorting allowed the high-dimensional phenotypic tracking of even modestly expanded B-lineage clones. To minimize PCR and sequencing errors, we used a high-fidelity PCR enzyme [Phusion High-Fidelity DNA Polymerase (NEB)] with an error rate > 50-fold lower than the typical Taq polymerases (41) and paired-end sequencing (27) that yielded hundreds to thousands of reads per cell allowing for correction of sequencing errors. With respect to sorting accuracy, the sorter identifies whether and to what extent markers are expressed on individual cells. However, fluorescence detection is subject to variance, which adds difficulty to assigning a cell to a population when that cell has characteristics that place it between two populations.

Our methodology integrates high-dimensional (13 fluorescence parameters plus forward-/side-scatter characteristics) phenotypic datasets into maps of multilevel intra- and interclonal diversity, which will support identification of disease-related phenotypes in autoimmune disorders, inflammation, infections, and malignancies beyond multiple myeloma.

Applying this methodology to B-lineage cells in multiple myeloma bone marrow, we also found the dominant multiple myeloma clones in up to 1.2% of sequenced cells that were phenotypically distinct from the majority. The existence of B-lineage cells other than plasma cells with the major clonotype and their physiological relevance has been discussed (19, 42–49). However, the phenotypic range of clonotypic cells could not be studied at the single-cell level in high-dimensional space. We show in multiple myeloma bone marrow that the phenotypes of clonotypic B-lineage cells can be heterogeneous, varying from differences in only 1 to 2 markers (e.g., expression of CD45 and absence of CD56 where the predominant clone is CD45⁻CD56⁺) to those similar to mature B cells (CD19⁺CD20⁺CD45⁺). Phenotypes of clonal cells were also heterogeneous, which would have given earlier FACS-sorting approaches with upfront-defined gating strategies limited chance of success. Given the diversity within the plasma cell compartment of our study samples and data in the literature (20, 21, 43, 44, 46, 50–54), low frequencies of different-phenotype B-lineage cells sharing the multiple myeloma rearrangement are not surprising. We would expect the variety of these phenotypes to increase with sample size. Whether these cells are multiple myeloma precursors, stem cells, or result from the phenotypic instability of multiple myeloma cannot be concluded from our data. However, the therapeutic success of anti-CD3 x anti-CD20 bispecific antibody-armed T cells (21) and chimeric antigen receptor T cells against CD19 in some patients with apparently CD19⁻ multiple myeloma (20) suggests that in selected cases, cells that are phenotypically B cells but harbor a plasma cell clonotype may be clinically relevant therapeutic targets. Our methodology could help to identify patients that might benefit from B-cell ablation.

Multiple myeloma patient 3 received treatment until two months before the bone marrow aspiration done for this study. Whether the lower phenotypic diversity of the cells belonging to the predominant multiple myeloma clone in this patient was due to the previous treatment or reflects the phenotypic distribution of the disease cannot be concluded from our data.

Future studies will determine whether clonotypic cells with B-cell immune phenotypes are more characteristic of mature B cells or of malignant plasma cells (13, 15, 16, 55). Our approach to identify these rare clonotypic cells requires sequencing, which in destroying the cells renders them inaccessible for downstream functional analyses.

Our approach also allowed us to discover other B-lineage clones that, although less abundant than the dominant species, were nonetheless amplified (between 1% and 9% of sequenced cells). Multiple B-lineage clones are evident in 2% to 6% of patients (34, 35) based on protein electrophoresis or immunofixation. These techniques have poor resolution, high detection limits (150–500 mg/L monoclonal protein; refs. 56–58), do not reach single-cell resolution, and require the monoclonal protein being present in the serum. With single-cell analysis, we found minor B-lineage clone expansion in all three of the patients analyzed. These minor clones showed phenotypes associated with mature B rather than plasma cell differentiation and seem to resemble physiologic processes as they also occur in healthy individuals.

Two patients shared the same CDR3 sequence as has been reported before in a patient with activation-induced cytidine deaminase deficiency (GenBank accession number EU789132; ref. 59). This and the fact that the predominant light chain CDR3 sequences from MM2 and MM3 have already been published in patients with multiple myeloma and cast nephropathy (60) suggests that there may be common antigen specificities in at least some patients with this disease.

The analysis of V gene mutation frequencies revealed varying and lower numbers of somatic mutations in cells belonging to 10 out of 12 of the less predominant clones. When comparing CDR3 nucleotide sequences within the same clones, we detected silent nucleotide exchanges, suggesting a convergent antigen-driven expansion of these clones as it also occurs independent form multiple myeloma in peripheral blood B cells from healthy individuals. Our methodology can give insights whether parallel expanding minor clones in a patient are likely to be part of the malignancy (aberrant phenotypes, nucleotide sequence identity including somatic mutations) and might be of potential therapeutic/prognostic relevance or follow disease-unrelated processes (normal B/plasma cell phenotypes, silent nucleotide exchanges).

In conclusion, we present a method that combines the accuracy of multiparameter FACS single-cell index-sorting with high throughput, high efficiency next-generation sequencing to link molecular information to immune phenotypes on the single-cell level. The method, applicable to the field of B-cell immunology, allows one to analyze the immune phenotypes that a given B- or plasma-cell clone can occupy, with the immunoglobulin sequence acting as a bar code to identify clonal progeny. We used this methodology for phenotypic tracking of B-lineage clones in multiple myeloma bone marrow and showed the expansion of molecularly defined predominant multiple myeloma clones to include (i) aberrant phenotype multiple myeloma cells and (ii) normal phenotype B-lineage cells. The parallel expansion of a variety of less predominant B-lineage clones includes (iii) aberrant phenotype, and (iv) normal phenotype B-lineage cells and in the majority of clones is the result of convergent B-lymphocyte expansion rather than malignant proliferation.

No potential conflicts of interest were disclosed.

Conception and design: L. Hansmann, M.M. Davis

Development of methodology: L. Hansmann, A. Han, M.M. Davis

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L. Hansmann, M. Liedtke

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L. Hansmann, L. Penter

Writing, review, and/or revision of the manuscript: L. Hansmann, M. Liedtke, M.M. Davis

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L. Hansmann, M.M. Davis

Study supervision: L. Hansmann, M.M. Davis

Other (figure design): L. Hansmann

We thank Francesco Vallania and Christopher Bolen for helpful discussions and support with bioinformatics, Hans-Peter Rahn at the Preparative Flowcytometry Facility at the Max Delbrück Center for Molecular Medicine for flow cytometry support, Kerstin Dietze for excellent technical assistance, and Alessandra Aquilanti for critically reading the manuscript. We thank Irene Panzer for the sharing of FACS reagents as well as Bernd Dörken and Jörg Westermann for their support. Single-cell sorting was performed in the Stanford Shared FACS Facility using NIH S10 Shared Instrument Grant (S10RR025518-01).

Leo Hansmann was supported by a research fellowship by the German Research Foundation (DFG, HA 6772/1-1) and is a Charité – Berlin Institute of Health (BIH) Clinician Scientist. Livius Penter is a BIH Junior Clinician Scientist. Major support for this work also came from the National Institutes of Health (U19 AI 057229 to M.M. Davis) and the Parker Institute for Cancer Immunotherapy.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.