The immune microenvironment shapes transcriptional and genetic heterogeneity in chronic lymphocytic leukemia

The compartmentalization of chronic lymphocytic leukemia (CLL) in peripheral blood (PB), bone marrow, and secondary lymphoid tissue has permitted the study of intrinsic and extrinsic drivers of tumorigenesis and progression. We previously showed that the lymph node (LN) tumor microenvironment (TME) activates B-cell receptor (BCR) and nuclear factor κ B signaling and induces the proliferation of CLL cells.¹ T–cell-derived CD40L and interleukin-4(IL), and possibly macrophage-derived B-cell activating factor (BAFF) and a proliferation-inducing ligand (APRIL), also promote tumor survival and signaling.^2-4 While receiving support from the TME, CLL cells modulate their surroundings by secreting chemokines, such as CCL3, CCL4, and CCL22, to attract T cells.⁵ These findings indicate complex, bidirectional crosstalk between CLL cells and their TME.

Beyond these microenvironmental dependencies, a diverse set of putative pathogenic mutations have been identified through large-scale sequencing efforts.^6,7 Heterogeneity in the genetic architecture of CLL has been closely linked to disease aggressiveness and its capacity to evolve under the selective pressures of treatment.^6,8 Clonal composition furthermore influences the CLL transcriptome in several ways. For example, biased expression of the mutated allele relative to the reference allele has been observed with mutations affecting DNA damage response genes such as ATM and TP53.⁹ Individual or cooccurrence of putative driver mutations can also result in more global changes in the transcriptome, dysregulating key oncogenic pathways.^10,11 

Herein, we integrate transcriptome profiling at the bulk and single-cell level with bulk whole-exome sequencing (WES) of paired PB and LN samples to understand the intersection between the TME and molecular events in CLL pathogenesis.

Details and resources are available in supplemental Methods.

Patient samples were obtained between 12 April 2005 and 26 August 2014 after written informed consent in accordance with the Declaration of Helsinki, applicable federal regulations, and requirements from the National Heart, Lung, and Blood Institute Institutional Review Board. The clinical study is registered at clinicaltrials.gov as NCT00923507.

PB mononuclear cells were isolated by density gradient centrifugation and cryopreserved. LN biopsies were either snap-frozen and stored at −80°C or mechanically disaggregated into single-cell suspensions (SCSs) and cryopreserved. CD19⁺ CLL cells were purified from PB mononuclear cells or LN SCSs by immunomagnetic positive selection with LS columns and stored as cell pellets at −80°C.

Total RNA libraries were prepared according to the Illumina TruSeq protocol and sequenced on a HiSeq 2000 (Illumina). Gene set enrichment analysis¹² was used to test for enrichment of a curated set of lymphocyte gene expression signatures (supplemental Table 2). CIBERSORT was used for deconvolution of bulk RNA-seq data from LN samples of patients with CLL and normal donors.¹ 

Single-cell libraries were created from LN SCSs with the Chromium Single Cell 3′ Reagent Kit V2. Libraries were sequenced on NovaSeq (Illumina) and processed by Cell Ranger Single Cell Software Suite 3.0.1.¹³ Downstream analyses were performed on Seurat v3.¹⁴ For RNA velocity, cell transitions were estimated based on quantification of unspliced and spliced mRNA.^15,16 

DNA was extracted from CD19⁺ CLL cells using DNeasy kits (Qiagen, Germantown, MD). Libraries were prepared using Agilent SureSelect Human All Exon kit and Nextera DNA Library Prep for Enrichment. Pooled libraries were sequenced on Illumina next-generation sequencers. Binary alignment map (BAM) files were generated from the Picard pipeline and then analyzed on the Firehose platform, including quality control, local realignment, mutation calling, small insertion and deletion identification, rearrangement detection, and coverage calculations. ABSOLUTE¹⁷ was used to estimate sample purity, ploidy, absolute somatic copy number, and the cancer cell fraction (CCF) of mutations. To distinguish subclonal populations, a Bayesian clustering procedure^18,19 was applied to cancer cell fractions (CCFs) across PB and LN samples from each patient.

Formalin-fixed, paraffin-embedded LN samples were stained with CD163. Images were captured at a magnification of 10× on an Olympus Bx41 microscope. The number of CD163⁺ cells was counted in 5 consecutive fields at 40× by a hematopathologist and then averaged to obtain the number of cells per high-power field.

Cryopreserved LN SCSs were thawed and subjected to CD19⁺ immunomagnetic depletion with LD columns. Cells were preincubated with LIVE/DEAD fixable Aqua Dead Cell Stain then stained with a cocktail of antibodies against CD3, CD4, CD8, CD45RO, CCR7, HLA-DR, CD19, and CD14. Fluorescence minus one controls were used to set gates for positive events. Cells were analyzed on LSRFortessa machine and FlowJo.

Continuous variables were compared between patient subgroups with the 2-sided Student t test. Correlations were computed using Pearson correlation. The Kaplan-Meier method was used to estimate progression-free survival, and differences between patient subgroups were compared using the log-rank test. Statistical analyses were performed in Prism 7 or R.

The bulk and scRNA-seq data have been deposited in the Gene Expression Omnibus database under accession numbers GSE161711 and GSE161610. The WES data have been deposited in the dbGaP database under dbGaP accession phs002297.v1.p1.

Prior studies of microarray-based gene expression profiling and in vivo deuterium labeling have revealed the distinct biology of CLL across anatomic compartments.^1,20 To deeply investigate the behavior of CLL in the TME, we performed total RNA sequencing (RNA-seq) of purified CD19⁺ tumor cells from paired PB and LN samples in 29 patients with treatment-naïve CLL or small lymphocytic lymphoma (Figure 1A,B; Table 1). Principal component analysis of bulk transcriptomes separated samples by tissue type (Figure 2A). We identified upregulation of 253 genes and downregulation of 32 genes in purified tumor cells from the LN compared with the PB (fold-change ≥2, false discovery rate [FDR] <0.05) (Figure 2B; supplemental Table 1).

To evaluate the impact of the TME on the cellular transcriptome, we applied gene set enrichment analysis using a curated list of gene expression signatures of signaling pathways, transcription factor regulation, and cellular processes such as proliferation and metabolism (Figure 2C; supplemental Table 2). Consistent with previously reported gene expression profiling,^1,21 we confirmed BCR, downstream NFκB, and NOTCH activation as enriched pathways in the LN. Crosstalk with the TME was supported by signatures of CD40 and cytokine signaling. Proliferation and metabolism signatures were also enriched in the LN. In contrast, PB tumor cells were characterized by 2 signatures: cellular quiescence and the transcription factor Krüppel-like factor 2, which maintains mature B cells in a resting state.²² 

Somatic hypermutation of the immunoglobulin heavy-chain variable region (IGHV) is associated with cellular energy manifesting as reduced signaling capacity upon in vitro stimulation of the BCR by anti-immunoglobulin–M and better prognosis.^23,24 To assess differences between mutated (M-CLL) and unmutated (U-CLL) CLL cases in the TME, we compared the expression of a BCR signature in LN samples. We found overexpression of this BCR signature in U-CLL compared with M-CLL (P = .03; unpaired t test) (Figure 2D). U-CLL cases also overexpressed NFκB, interferon, and T-cell cytokine signatures (P < .05). No signatures were overexpressed in M-CLL cases.

To further dissect the transcriptional heterogeneity of CLL cells within the TME, we performed a high-throughput, droplet-based single-cell RNA-seq of 5 LN samples. We sequenced a median of 2946 (range, 2605-3377) cells per sample and detected a median of 1171 (range, 916-1374) genes per cell. A total of 15 107 cells were analyzed by Seurat v3¹⁴ using anchor correspondences across samples to distinguish shared and nonoverlapping cell states (Figure 3A). Three cell lineages (CLL, T cells, and monocytes) were identified (Figure 3B,E). CLL cells were comprised of 3 major subpopulations defined by quiescent, activated, and proliferating cell states (Figure 3B,E). Compared with the major subpopulation of quiescent CLL cells, activated CLL cells, accounting for 2.2% to 4.3% of tumor cells, overexpressed BCR, NOTCH, proliferation, and metabolism signatures (Figure 3F). Proliferating CLL cells, identified by G2/M and S phase markers, comprised 0.4% to 1.0% of tumor cells (supplemental Figure 1A).

To infer directional dynamics of the 3 CLL cell states, we estimated RNA velocity based on the relative abundance of unspliced pre-mRNA and spliced mRNA in single cells.^15,16 For a given gene, unspliced mRNA corresponds to the future state of spliced mRNA. When the balance of unspliced and spliced mRNAs of many genes are analyzed in each single-cell transcriptome, the RNA velocity of these cells can be estimated, which can, in aggregate, reveal the movement from one cell state to the next. In our dataset, the dominant structure comprised of proliferating cells as the root, giving rise to activated cells before entering the final state as quiescent cells (Figure 3C). This program, consistently observed across 3 different types of RNA velocity modeling (supplemental Figure 1B), implies that CLL cell fate in the LN is unidirectional and that most quiescent cells are not reactivated by TME stimuli.

We reasoned the bulk transcriptional profiles of PB and LN samples might be related to the tissue-specific distribution of CLL cell states. Therefore, scRNA-seq was also performed in 5 matched PB samples. A median of 5754 (range, 5071-6248) cells per sample and 29 008 cells in total were projected onto the LN scRNA-seq data to quantify the proportion of tumor cells in each cell state (Figure 3D; supplemental Figure 2A). Proliferating CLL cells were nearly absent in PB, and activated cells were significantly less abundant in PB than LN (Figure 3G). Furthermore, the activated population in PB underexpressed cell cycle genes, CCND2 and DUSP2, and chemokine CCL4 compared with their counterpart in LN (P < .05 for all comparisons) (supplemental Figure 2B).

Using single-cell expression profiles, we derived an activated CLL signature of the 10 most overexpressed genes by mean fold-change (HSP90AB1, ENO1, TUBB, RAN, PRDX1, LDHA, NME1, DNPH1, HSPD1, PKM) among activated cells. We assessed the expression of this signature in bulk RNA-seq of unselected LN samples from 34 patients. These 34 patients included 29 patients who contributed paired LN and PB samples to the tumor transcriptome analysis above and 5 additional patients with small lymphocytic lymphoma (Figure 1A). The activated CLL signature was overexpressed in U-CLL compared with M-CLL (P = .0033), in younger patients <60 years of age compared with older patients (P = .014), and in CLL compared with small lymphocytic lymphoma (P = .049, all unpaired t tests) (Figure 4A). No differences were observed among cytogenetic subgroups (Figure 4B).

Patients were followed for a median of 62.5 (range, 0.5-155.7) months after their LN biopsy. Consistent with more aggressive disease, higher expression of the activated CLL signature was associated with inferior treatment-free survival (TFS) (hazard ratio [HR], 3.9; 95% confidence interval [CI], 1.7-8.9; P < .0001) (Figure 4C). Median TFS from the time of diagnosis to first treatment or death was 29.0 months among patients with high signature expression and 124.6 months among those with low signature expression. In comparison, we analyzed TFS based on IGHV mutational status, a well-established prognostic marker in CLL.^25,26 As expected, TFS of U-CLL was inferior to M-CLL, with a median TFS of 29.2 months vs 137 months (HR, 4.0; 95% CI, 1.9-8.5; P = .0002) (Figure 4D). Within each IGHV subgroup, the activated CLL signature retained independent prognostic information. Median TFS in patients with M-CLL was 76.5 months with high signature expression vs 186.5 months with low signature expression (HR, 3.8; 95% CI, 0.7-20.7; P = .040) (Figure 4E), and in patients with U-CLL was 16.1 months with high signature expression vs 37.4 months with low signature expression (HR, 2.5; 95% CI, 1.0-6.2; P = .025) (Figure 4F).

Although cellular interactions within the TME are thought to be critical in the pathogenesis of CLL, understanding the contribution of nonmalignant immune cells has been largely based upon in vitro and murine models. Whereas most studies have focused on specific cell types of interest, we sought to comprehensively assess the immune cell repertoire within the LN TME. The proportion of various cell types can be inferred by deconvolving bulk gene expression profiles.^27,28 CIBERSORT^27,28 was used to construct a custom gene signature matrix from gene expression profiles of purified CLL cells from LN samples and 20 non-B immune cell types. This signature matrix was then applied to deconvolve bulk RNA-seq data of the 34 unselected CLL LN samples. In parallel, we performed bulk RNA-seq of 4 normal LN samples and deconvolved these data using a validated gene signature matrix of 22 immune cell types, including naïve and memory B cells.²⁸ The median estimated proportion of tumor cells was 77.9% (range, 29.6% to 90.7%). Only 1 sample had <50% tumor content. B cells comprised a median of 34.5% (range, 28.5% to 42.4%) of normal LNs and were predominantly memory B cells. After normalization to tumor or B-cell content, we compared the proportions of the 20 non-B immune cell types between CLL and normal LNs (Figure 5A; supplemental Figure 3). Follicular helper T cells, which provide T-cell help to B cells within germinal centers,²⁹ were expanded while naïve CD4⁺ T cells and uncommitted macrophages were decreased in CLL LNs compared with normal LNs.

We hypothesized that a subpopulation of tumor cells, specifically that of the activated phenotype, could be interacting with immune cells of the TME. Expression of the activated CLL signature among 34 CLL LNs was correlated against the proportion of 20 non-B immune cell types inferred by deconvolution. The activated CLL signature was significantly correlated with M2 macrophages (Pearson R = 0.44; P = .0087) (Figure 5B) and activated CD4⁺ memory T cells (Pearson R = 0.53; P = .0013) (Figure 5D). We corroborated these in silico findings with immunohistochemical staining of CD163 for M2 macrophages (Figure 5C) and flow cytometric analysis of activated HLA-DR⁺ activated CD4⁺ effector and central memory T cells in a subset of CLL LNs (Figure 5E; supplemental Figure 4). Given the association between the activated CLL signature and IGHV mutation status, we compared the TME immune cell composition between U-CLL and M-CLL. The proportion of the 20 non-B immune cell types was not significantly different based on IGHV mutation status. These data indicate that the TME is abnormal in CLL LNs compared with normal LNs and that an activated subset of CLL cells shapes the immune cell milieu in which they reside.

Genetic lesions have been identified exclusively in either the PB or LN in some patients with CLL.³⁰ We sought to map genotype–phenotype relationships that could explain the compartmentalization of subclones. WES was performed on CD19⁺ cells of paired PB and LN samples and matched germline DNA from 14 patients. A median of 27 (range, 11-69) somatic single nucleotide variants (sSNVs) and 3 (0-10) insertions and deletions (sIndels) were detected per exome. There were 28 sSNVs and sIndels shared by ≥2 patients (Figure 6A). In addition, all but 1 patient had copy number alterations, most commonly del(11q) and del(13q).

Cancer cell fractions (CCFs) of sSNVs, sIndels, and copy number alterations were inferred from variant allele frequencies using ABSOLUTE,¹⁷ then clustered over the 2 anatomic compartments for each patient. Genetic compartmentalization (ΔCCF >0.25; FDR <0.1) was observed in 7 patients (50%), of whom 6 demonstrated subclonal expansion in the LN (Figure 6B). In the remaining 7 patients, no genetic compartmentalization was observed (supplemental Figure 5). To explore factors that contribute to the enrichment of subclones in the TME, we grouped patients based on the presence (shifted group) or absence (stable group) of an expanded subclone in LN. Patient CLL-C46 showed subclonal expansion in PB and was excluded.

Activation-induced cytidine deaminase (AID) is required for somatic hypermutation and class switch recombination but can also give rise to oncogenic mutations.³¹ In the TCL1 model, AID action has been implicated in clonal evolution.³² Therefore, we assessed AICDA expression in bulk RNA-seq data of purified CLL cells from LN and PB. As expected, LN-resident CLL cells overexpressed AICDA compared with their circulating counterpart (P < .0001) (supplemental Figure 6A). Among LN samples, AICDA expression was significantly higher in U-CLL than in M-CLL (P = .0004) (supplemental Figure 6B). AICDA expression was weakly correlated with the activated CLL signature in bulk RNA-seq (supplemental Figure 6C). However, assessment of AICDA in proliferating, activated, and quiescent single cells could not be performed because AICDA transcripts were detected in only 0.84% of CLL cells in scRNA-seq. Notably, there was no difference in AICDA expression between the shifted and stable groups (supplemental Figure 6C).

Differential gene expression analysis between PB and LN was performed separately in the shifted and stable group. Most differentially expressed genes between PB and LN were shared by both groups of patients. However, each group preferentially overexpressed a small subset of genes in LN relative to PB (log₂FC >0.5; FDR <0.05) (Figure 6C). We observed upregulation of cell cycle genes (eg, E2F2, CDC25A, ESPL1) in the shifted group (supplemental Table 3). In the stable group, CD69, a lymphocyte activation marker, and CD83, involved in antigen presentation,³³ were upregulated, raising the question of whether immune-mediated control could restrict clonal outgrowth.

To explore this possibility, we used a tumor immune signature comprised of 18 T–cell-associated inflammatory genes, which was developed to predict the clinical outcome with checkpoint blockade in multiple cancer types.³⁴ Unsupervised hierarchal clustering of these signature genes in LN bulk RNA-seq data revealed increased expression in the stable group compared with the shifted group (Figure 6D). Irrespective of the source of mutagenesis, our findings raise the possibility of T–cell-mediated immune surveillance as a factor restricting clonal evolution in the LN.

The development and progression of CLL are not isolated to intrinsic properties of leukemia cells but rather involve complex interactions with the TME. In vitro experiments and murine models have identified candidate cellular and soluble components in the local milieu that may be coopted by tumors to support leukemogenesis. In this study, we integrated bulk and single-cell transcriptome and bulk exome sequencing of paired PB and LN samples to define the in situ interactions between CLL and the TME and map clonal trajectories between anatomic compartments.

Our bulk RNA-seq data provided a comprehensive view of signaling pathways and cellular processes that are modulated in leukemic cells by the TME. These findings indicate that multiple support signals from the TME collectively promote growth and survival. Our sample size was sufficient in informing the differences between M-CLL and U-CLL in the TME. Compared with M-CLL, U-CLL demonstrated not only enhanced BCR signaling, consistent with prior studies,²³ but also increased engagement with immune cells of the TME, as evidenced by the overexpression of interferon and T-cell cytokine signatures. While both M-CLL and U-CLL exhibit autonomous BCR signaling,³⁵ our data demonstrate that extrinsic factors in the TME augment BCR activation above the basal level established by autonomous signaling in U-CLL more so than in M-CLL.

The gene expression profile provided by bulk RNA-seq is averaged across many cells, each of which has an individual transcriptome. By performing scRNA-seq in a subset of paired PB and LN samples, we found that bulk transcriptome data are, in fact, the aggregate of distinct CLL subpopulations. The transcriptional program associated with LN-resident CLL, that is, increased signaling and proliferation compared with its counterpart in PB, is derived from minor subpopulations comprising <10% of all tumor cells in the TME. By analyzing the RNA velocity of single cells, we found that CLL cells start in the proliferating state, move to the activated state, and end in the quiescent state. The sequential progression of these cell states suggests that most quiescent CLL cells residing in the LN are never reactivated by TME stimuli.

In contrast to the insights gained about tumor heterogeneity, we were only able to identify 3 nontumor immune cell clusters by scRNA-seq of LN samples. Several factors may have contributed to this result. First, adherent populations such as dendritic cells and macrophages are lost during the mechanical dissociation of LN tissue into single cells. Second, the relatively low number of transcripts per cell, which is inherent to scRNA-seq, limits the resolution of closely related immune cell subsets. To overcome these limitations and to analyze more samples in a cost-effective manner, we used in silico transcriptome deconvolution to profile accessory cell types in the LN TME of patients with CLL compared with normal LN. Our data indicate the activated subpopulation shapes the local immune cell composition, specifically of M2 macrophages and activated CD4⁺ memory T cells, independently confirming conclusions derived from model systems.^36-38 The proportion of activated CLL cells and their capacity to recruit a tumor-supportive TME may underlie the shorter TFS observed in patients with higher expression of the activated CLL signature. Taken together, an identifiable minor tumor subpopulation underlies CLL pathogenesis and drives clinical outcomes.

A prior study detected PB- and LN-specific genetic lesions in 6 of 9 patients,³⁰ but this interpretation overlooked possible differences in the size of subclones common to both PB and LN. In contrast, we included shared and specific mutations to infer and compare the clonal architecture between anatomic compartments. In one-half of patients, the CCF of identical subclones differed by >25% between PB and LN. Although certain mutations could, in theory, confer relative clonal fitness in a given environment, perhaps limited by sample size, we did not identify specific mutations that were preferentially distributed in either PB or LN. Together with the RNA velocity of CLL single cells, the spatial imbalance in genetic composition challenges the lifecycle model of CLL, in which tumor cells transit from a proliferative state in lymphoid tissue to a resting state in the blood and back.³⁹ 

Our data support differential fitness of subclones within the TME driving genetic compartmentalization. The immune TME is capable of constraining clonal expansion and maintaining tumors in a state of equilibrium.⁴⁰ In previously treated CLL, nodal relapse without leukemic progression implicates the role of protective niches in clonal evolution.⁴¹ Consistent with a permissive TME, subclonal expansion was apparent in the LN in all but 1 case with genetic compartmentalization. Bulk RNA-seq data from the same samples suggested that clonal outgrowth is impeded by an activated immune TME. When immune surveillance is inactivated, subclones with a competitive advantage may expand in response to support signals provided by the TME.

Experimental models have implicated AID in promoting clonal evolution.³² We observed increased AICDA expression in LN compared with PB and in U-CLL compared with M-CLL. However, there was no difference in AICDA expression between cases with or without clonal evolution in situ. Irrespective of the source of mutagenesis, our findings indicate T–cell-mediated immune surveillance as a factor that may restrict clonal evolution in the LN. Cumulative evidence has pointed to the existence of CLL-specific T cells and treatment, such as ibrutinib or lenalidomide, that could restore their cytotoxic activity.^42-44 Therefore, strategies to strengthen host immune surveillance in CLL could represent one avenue to impede the evolution to more aggressive clones. Collectively, our findings underscore the value of an integrated multiomic approach to the evaluation of the LN TME in CLL.

The authors thank the patients, the National Heart, Lung, and Blood Institute (NHLBI) Bioinformatics and Computational Biology Core, the NHLBI Sequencing and Genomics Core, and J. Joseph Melenhorst for his expertise in flow cytometric analysis of T-cell subsets.

This research was supported by the Intramural Research Program of the NHLBI, National Institutes of Health.

Contribution: C.S. and A.W. conceptualized the study; C.S., Y.-C.C., A.M.Z., M.J.B., S.P., D.L., D.R., and S.G. provided formal analysis; C.S., M.J.B., N.S.S., and T.D.-H. undertook investigation; C.S. and A.W. provided resources; C.S. curated data; C.S. and A.W. wrote the original manuscript draft; C.S., Y.-C.C., and A.M.Z. performed visualization; G.G., M.P., C.J.W., and A.W. supervised the project; C.J.W. and A.W. acquired funding; and all authors reviewed and edited the manuscript.

Conflict-of-interest disclosure: C.S. receives research support from Genmab. M.J.B. is an employee at AstraZeneca. S.H.G. is a consultant for Novalgen Ltd and holds patents related to ROR1 therapies; has received honoraria, conference support, and speakers fees from AstraZeneca, AbbVie, and Janssen. N.S.S. is a consultant for and has received speaker fees from AbbVie, Pharmacyclics LLC, an AbbVie company, and Janssen; and is a member of the advisory board for TG Therapeutics, Innocare, BeiGene, KyowaKirin, ADC therapeutics, and Kite. G.G. receives research funds from IBM and Pharmacyclics; and is an inventor on patent applications related to MSMuTect, MSMutSig, MSIDetect, POLYSOLVER, and SignatureAnalyzer-GPU; is a founder, consultant, and holds privately held equity in Scorpion Therapeutics. C.J.W. holds equity in BioNTech, Inc and receives research funding from Pharmacyclics. A.W. receives research support from Pharmacyclics LLC, an AbbVie Company, Acerta LLC, a member of the AstraZeneca Group, Merck, Verastem, Nurix, and Genmab. The remaining authors declare no competing financial interests.

Correspondence: Clare Sun, Hematology Branch, NHLBI, NIH, Building 10, CRC 3-5132, 10 Center Dr, Bethesda, MD 20892-1202; e-mail: clare.sun@nih.gov; and Adrian Wiestner, Hematology Branch, NHLBI, NIH, Building 10, CRC 3-5140, 10 Center Dr, Bethesda, MD 20892-1202; e-mail: wiestnea@nhlbi.nih.gov.