Single-cell CRISPR screening characterizes transcriptional deregulation in T-cell acute lymphoblastic leukemia

Abstract

T-cell acute lymphoblastic leukemia (T-ALL) is an aggressive type of leukemia caused by accumulation of multiple genetic alterations in T-cell progenitors. However, for many genes it remains unknown how their mutations contribute to disease development. We therefore performed two single-cell CRISPR screens in primary pro-T cells ex vivo to study the transcriptional impact of loss-of-function alterations in T-ALL and correlate this with effects on cell fitness. The various perturbations were clustered based on their effects on E2F/MYC or STAT/NOTCH signatures, which play a defining role in driving T-cell proliferation. Many of the perturbations resulted in positive effects on the STAT and NOTCH signatures and were predicted to behave as haploinsufficient tumor suppressors in T-ALL. Additionally, Spi1 was identified as an essential gene for pro-T-cell survival, associated with deregulation of the MYC signature and epigenetic consequences. In contrast, Bcl11b was identified as a strong tumor suppressor gene in immature T lymphocytes, associated with deregulation of NF-kB and JAK/STAT signaling. We found a correlation between BCL11B expression level and JAK/STAT pathway mutations in T-ALL patients and demonstrated oncogenic cooperation between Bcl11b inactivation and JAK3 hyperactivation in pro-T cells. Altogether, these single-cell CRISPR screens in pro-T cells provide fundamental insights into the mechanisms of transcriptional deregulation caused by genetic alterations in T-ALL.

Introduction

T-cell acute lymphoblastic leukemia (T-ALL) is an aggressive form of leukemia characterized by high white blood cell counts and infiltration of immature T cells into the bone marrow and other tissues. Sequencing studies of patient cohorts have identified many genetic alterations in T-ALL patients, often in major signaling pathways regulating T-cell development such as NOTCH, JAK/STAT, PI3K/AKT, RAS and the T-cell receptor signaling pathway.^1-4 While frequent oncogenic events such as NOTCH1 mutations have been extensively studied, the contribution of other mutations, for example in epigenetic regulators such as CTCF, remains poorly understood.²

Many of the current in vitro models for T-ALL have their limitations, since immortalized leukemic cell lines intrinsically have a complex genetic background with numerous mutations and structural abnormalities. In recent years, a pro-T-cell model was developed to culture murine primary T-cell precursors ex vivo.⁵ This model has a wild-type genetic background, facilitating the study of gene functions and oncogenic cooperation in T-ALL in a clean way. Pro-T cells are generated by isolation of murine lineage-negative bone marrow cells and culturing these primary cells ex vivo in the presence of murine DLL4 (stimulating NOTCH), murine IL7 (stimulating JAK/STAT) and murine SCF (stimulating c-KIT). This drives the differentiation of the primary cells into immature double-negative T lymphocytes, which were shown to represent a suitable model for studying T-cell development and transformation.^6-8 We have previously shown that recurrent mutations in T-ALL can drive cytokine-independent proliferation of pro-T cells, demonstrating that this model effectively simulates leukemic transformation ex vivo.^6-10

Because many of the recurrent genetic alterations in T-ALL patients remain largely uncharacterized, we used single-cell CRISPR (clustered regularly interspaced short palindromic repeats) screening to study a number of potential tumor suppressor genes in pro-T cells. This technique combines CRISPR/Cas9 screening with single-cell RNA sequencing and allows the interrogation of multiple genetic alterations while studying the full transcriptome associated with each perturbation.^11-14 Furthermore, enrichment or depletion of the guide RNA (gRNA) can be assessed to evaluate the impact of each perturbation on cell fitness. Specific measures are required for detection of the gRNA in a single cell, either through a unique barcode^13-15 or by directly sequencing the gRNA using a modified construct with a poly-A tail¹² or capture sequence.¹⁶ Single-cell CRISPR screening has been applied previously in T cells by Zhou et al.,¹⁷ who studied different transcription factors that drive the T-cell developmental trajectory. In general, single-cell sequencing is increasingly used to study the heterogeneity of leukemia at single-cell resolution.^18-20

Alongside a number of major T-cell regulators that are known to play a role in T-ALL development, we targeted multiple genes that are recurrently mutated in patients but for which their exact function remains to be elucidated. In a first experiment, we perturbed 17 known transcriptional regulators, to enable optimalization of the methodology and identification of the driving gene signatures. We subsequently designed a larger library targeting 42 recurrently mutated genes, with the aim of characterizing their transcriptional profiles, providing fundamental insight into the function of these genes and the affected downstream signaling pathways.

Methods

Experimental design

Single-cell CRISPR screening was performed on ex vivo-cultured pro-T lymphocytes to study recurrent mutations in T-ALL. Custom gRNA libraries were designed and retrovirally transduced in pro-T cells, with subsequent fluorescence-activated cell sorting (FACS) and cell harvesting at multiple timepoints. Samples were loaded onto the 10X Genomics Chromium for single-cell encapsulation and mRNA and gRNA extraction. Transcriptome data were analyzed to obtain gene expression signatures per perturbation and gRNA counts were used to infer proliferative effects.

Cell culture and transduction

Pro-T-cell cultures were established as described previously⁵ from C57BL/6J (Charles River Laboratories, Saint-Germain-Nuelles, France) Cas9 knock-in transgenic mice (024858, Jackson Laboratories, Bar Harbor, ME, USA). Lineage-negative hematopoietic stem cells were purified from the bone marrow of C57BL/6 mice using an EasySep Mouse Hematopoietic Progenitor Cell Isolation Kit (StemCell Technologies). The use of mice for the purpose of establishing pro-T-cell cultures was approved and supervised by the KU Leuven ethical committee (ECD 030/2023, approval date March 16, 2023). Pro-T cells were maintained in non-tissue culture-treated plates pre-coated with anti-human IgG1-Fc monoclonal antibody (Abcam ab1927). Pre-coated wells were washed twice and subsequently coated with 2 μg/mL murine DLL4-Fc in RPMI-1640. Pro-T cells were seeded at 500,000 cells/mL in medium supplemented with 20% fetal bovine serum with murine IL7 (20 ng/mL, PeproTech), murine SCF (20 ng/mL, PeproTech) and primocin (100 μg/mL, Invivogen). Cells were transduced by retroviral spin-fection at 2,500 rpm for 90 minutes at 30°C, with 0.5 μL/mL polybrene (10 mg/mL, Sigma Aldrich TR-1003-G) per 10⁶ cells. Next, cells were incubated at 37°C for 2 hours, washed with phosphate-buffered saline and seeded into freshly coated wells. For proliferation assays using gRNA or short hairpin (sh)RNA, cells were seeded at 500,000 cells/mL and cell number, viability and fluorescent signal were monitored using a MACSQuant VYB (Miltenyi) flow cytometer. For drug treatments, 50,000 cells were seeded in 100 μL RPMI-1640 medium with 20% fetal bovine serum with murine IL7, murine SCF and primocin in 96-well plates. The D300e digital dispenser (TECAN) was used for dispensing specific concentrations of ruxolitinib phosphate (MedChemExpress) in triplicate. Cell viability was evaluated after 60 hours using the ATPlite luminescence assay (PerkinElmer) on a Victor X4 plate reader (PerkinElmer).

Single-cell CRISPR screening

Pro-T cells with constitutive Cas9 expression were transduced with the CRISPR library with an efficiency of around 10%. Transduced cells were sorted 48 hours after transduction using an Influx cell sorter (BD Biosciences) or FACSAria Fusion (BD Biosciences). Sorted cells were cultured for a total of 14 days, with harvesting at days 3, 7 and 14. At each timepoint, cell suspensions of 1,000 cells/L in phosphate-buffered saline + 0.04% bovine serum albumin were loaded onto the Chromium controller (10X Genomics) for single-cell encapsulation. Sequencing libraries were prepared according to the manufacturer’s directions using Chromium Next GEM single-cell 3’ reagents with feature barcoding technology. Libraries were analyzed with a High Sensitivity DNA kit on a Bioanalyzer (Agilent Technologies) and sequenced on a NovaSeq 6000 instrument (Illumina).

Statistical analysis and data visualization

Graphpad Prism (version 9.2.0) was used for statistical analyses. Data were displayed as mean ± standard deviation. Comparisons between two groups were performed by unpaired two-tailed Student t tests. One-way analysis of variance (ANOVA) was used for comparisons of multiple groups in proliferation assays, applying the Tukey correction for multiple testing. Prior to ANOVA, a Shapiro-Wilk normality test was performed to test for Gaussian distribution of residuals and Brown-Forsythe test for homogeneity of variance. The statistical significance of P values is indicated by asterisks (****P<0.0001, ***P<0.001, **P<0.01, *P<0.05) or NS (not significant). Significance of overlap between two lists was calculated via the hypergeometric probability distribution and visualized in Venn diagrams. Box plots show a line at the median with the box extending from the first (Q1) to third (Q3) quartiles, with whiskers ranging from minimum to maximum values. FlowJo (version 10.8.1) was used for analysis of flow cytometry data. Custom R (version 4.1.1) scripts were used to analyze single-cell and bulk RNA-sequencing data, as well as for data visualization, as further elaborated on in the Online Supplementary Methods. Webgestalt²¹ and the Broad GSEA software (version 4.2.1)^22,23 were used for analysis of differentially expressed genes, and i-cisTarget²⁴ for prediction of cis-regulatory motifs.

Results

Transcriptional phenotyping of perturbations in pro-T cells by single-cell CRISPR screening

We first performed a CRISPR screen using a gRNA library targeting 17 important T-cell regulators with a known role in T-ALL (Online Supplementary Table S1). Five gRNA were designed per gene, together with three non-expressed genes as negative controls (NT1-3), resulting in a library of 96 gRNA (Online Supplementary Table S1). The 10X Genomics Chromium 3’ RNA-sequencing chemistry with feature barcoding technology was used to detect gRNA via their capture sequence (Online Supplementary Figure S1A).¹⁶ Different capture sequence configurations were evaluated experimentally first in Ba/F3 cells and subsequently in the pro-T-cell model (Online Supplementary Figure S1B).

Pro-T cells were transduced with the gRNA library with a low multiplicity of infection (<10%) to avoid introducing multiple gRNA in a single cell. The transduced cells were subsequently selected by FACS (Online Supplementary Figure S1A, C) and samples were harvested at days 3, 7 and 14 after transduction (Figure 1A). After filtering bad quality cells and doublets, the dataset contained 4,773, 3,606 and 9,672 cells, for each respective timepoint, with 45% to 54% of cells containing a single gRNA (Figure 1B). At day 3, a median of 112 good-quality cells (interquartile range [IQR]: 97-133) per target gene and a median of 21 cells (IQR: 16-27) per gRNA were sequenced (Figure 1C). At this timepoint, we observed an almost equal distribution of gRNA abundance across the cell population (Figure 1D, Online Supplementary Figure S1D). Since the cells were sequenced at multiple timepoints, the impact of each perturbation on cell proliferation could be studied, thereby identifying essential genes and tumor suppressors. Cells with perturbation of Il7r, Notch1, Myc, Kit and Spi1 had a significant proliferative disadvantage, in agreement with the fact that pro-T cells are dependent on JAK/STAT, NOTCH and KIT stimulation.⁵ In contrast, cells with Pten inactivation had a strong growth advantage (Figure 1D, E), which was expected since PTEN is a well-described tumor suppressor.²⁵ Cells with Pten gRNA had become the dominant clone by day 14, covering over 51% of the cell population and outcompeting clones with other perturbations. Comparison of the trends of all gRNA over time with the NT cells also identified Fos, Jun, Ptprc, E2f1 and Runx1 as mild tumor suppressor genes (Figure 1E).

Perturbation of major transcriptional regulators identifies essential gene expression signatures

To study transcriptome profiles at the single-cell level, the cDNA libraries were sequenced with 108,400, 160,011 and 59,678 mean reads per cell for days 3, 7 and 14, respectively. This resulted in a median number of detected transcripts per cell of 4,285 (IQR: 3,406-5,264), 5,106 (IQR: 4,186-5,889) and 3,649 (IQR: 2,948-4,390), respectively (Figure 2A). We observed efficient CRISPR-mediated downregulation for the majority of target genes (Figure 2B, Online Supplementary Figure S2). However, downregulation of gene expression is not a perfect measure for knockdown efficiency as CRISPR editing could result in intact mRNA generating a non-functional protein. Indeed, for some genes, such as Myc and Il7r, we did not observe downregulation of the gene itself, but we did detect changes in expression of known downstream target genes (Online Supplementary Figures S1E and S2). Cis-regulatory motif analysis of the significantly differentially expressed genes also demonstrated enrichment of these transcription factor binding motifs (Figure 2C). In summary, combined analysis of the expression of the gene targeted by the gRNA, the expression of downstream targets and a regulatory motif analysis provides a robust analysis to evaluate perturbation efficiency in silico. Cell cycle phase was estimated based on the expression of cell cycle marker genes. A strong correlation was observed between proliferation based on cell number (Figure 1E) and predicted phase distribution based on marker gene expression (Figure 2D). In cells with Myc, Kit, Notch1 and Spi1 gRNA, a decreased fraction of cells resided in the G2M/S phase, while the number of cells in the stationary G1 phase was elevated. In contrast, cells with Pten, E2f1, Ptprc, Cdkn2a, Jun and Runx1 inactivation showed a higher proportion of cells in the G2M/S phases, indicative of more active cycling (Figures 1E and 2D).

Based on the single-cell RNA-sequencing data, we aimed to identify core gene expression signatures driving the response of pro-T cells to CRISPR perturbation. Since T lymphocytes are highly dependent on STAT and NOTCH signaling, with MYC as a downstream effector, we focused on these pathways in particular. Enrichment analyses were performed on the ranked lists of differentially expressed genes of the pseudo-bulk transcriptomes, using a selection of JAK/STAT, NOTCH and MYC gene sets.^7,26,27 Significant negative enrichment of multiple relevant gene sets was observed in cells with the respective perturbations. The signatures described by Bornschein et al.⁷ ranked among the most significantly enriched gene sets and were presumed to describe the behavior of pro-T cells most accurately (Table 1, Figure 2E).

Spi1 is an essential gene linked to MYC transcriptional activity

PU.1, encoded by Spi1, is an essential transcriptional regulator responsible for regulation of both myeloid and lymphoid development, in a time- and dose-dependent manner. In the T-cell lineage, Spi1 is typically only expressed in early progenitors and is responsible for delaying commitment and maintaining a multipotent state while allowing strong proliferation of immature T-cell precursors. Upon commitment, Spi1 is silenced, which is associated with upregulation of specific T-lineage genes and further progression through T-cell development.^28,29 Our CRISPR data identified Spi1 as an essential gene for pro-T cell proliferation (Figure 1E). Independent proliferation assays were performed to confirm the negative impact of Spi1 inactivation on cell fitness (Figure 3A, B), while contrariwise, Spi1 overexpression had no impact on proliferation (Figure 3C). To complement the single-cell RNA-sequencing data, bulk RNA sequencing was performed on cells with Spi1 inactivation. We found 3,097 significantly differentially expressed genes upon Spi1 loss, with a significant overlap between the single-cell and bulk datasets (Figure 3D, Online Supplementary Figure S3A, Online Supplementary Table S2).

Spi1 inactivation was associated with a significant negative enrichment of MYC target genes in both the bulk data as well the single-cell dataset (Figure 3E, F). Furthermore, cell cycle regulators, glycolysis and oxidative phosphorylation were downregulated, as well as E2F target genes (Figure 3E). Single-Cell regulatory Network Inference and Clustering (SCENIC)³⁰ was used to infer regulons associated with Spi1 perturbation based on the bulk data. This analysis also showed downregulation of Myc-related regulons, as well as downregulation of the E2f and Cebpb regulons, the latter being a known SPI1 binding partner (Figure 3G, Online Supplementary Figure S3B, C). Since Myc and its downstream targets are crucial for pro-T-cell survival and since the E2F transcription factor family is heavily involved in cell cycle regulation,³¹ downregulation of these signatures is in agreement with the observed negative impact of Spi1 inactivation on cell proliferation. In contrast to the negative effects on Myc, regulons associated with epigenetic regulation as well as regulons associated with Ets1 and other ETS factors were upregulated, indicating an impact on the epigenetic landscape and the ETS transcription factor family upon Spi1 loss (Figure 3G, Online Supplementary Figure S3D, E).

Table 1.Gene set enrichment analysis for cells with Il7r and Stat5b gRNA, Notch1 gRNA and Myc gRNA.

Figure 1.Pilot experiment of single-cell CRISPR screening in ex vivo pro-T cells. (A) Experimental setup of the single-cell CRISPR screens. Pro-T cells were generated by isolation of lineage-negative hematopoietic stem and progenitor cells from C57BL/6J CAS9 knock-in mice which were cultured ex vivo in the presence of murine DLL4, SCF and IL7. Custom gRNA libraries were designed and cloned in a retroviral construct with a capture sequence. After transduction and selection by fluorescence activated cell sorting, cells were cultured for a total of 14 days, with cell harvesting at days 3, 7 and 14 for single-cell library preparation of the mRNA and gRNA. (B) Pie charts show percentages of cells with zero, one or multiple gRNA, with indication of the total number of cells passing filtering parameters at each timepoint. (C) Number of cells per gene and per gRNA. The gRNA library contained five gRNA per target gene. Box plots show the number of sequenced cells passing filter criteria per gene (upper panel) and per gRNA (lower panel). Plots show a line at the median and the box extends from the first quartile (Q1) to the third quartile (Q3), with whiskers ranging from minimum to maximum values. (D) Distribution of cells per CRISPR target gene at days 3, 7 and 14. Each slice of the donut represents the fraction of cells with a specific perturbation. Over time, changes in the ratio are observed due to differences in proliferation rate in cells with different perturbations. (E) Evolution of the fraction of cells per perturbation over time across days 3, 7 and 14. The Y-axis shows data on a logarithmic scale. Non-targeting controls (NT1, NT2, and NT3) are indicated in gray font. Perturbations are ordered based on the number of cells at day 14. gRNA: guide RNA; Lin-: lineage-negative; HSPC: hematopoietic stem and progenitor cells; mIL7: murine interleukin-7; mSCF: murine stem cell factor; mDLL4: murine DLL4; FACS: fluorescence activated cell sorting.

Figure 2.Single-cell RNA-sequencing data demonstrated downregulation of target genes, downstream pathways and predicted cis-regulatory motifs. (A) Violin plots displaying the number of detected transcripts per cell for each timepoint. (B) Normalized expression levels of the target genes in cells with the respective gRNA. To study gRNA efficiency, all cells from the three timepoints were aggregated per perturbation to create ridge plots of the normalized expression level of each target gene. For which transcript the expression level is displayed is indicated at the top of the graph, with the purple ridge highlighting cells with the corresponding perturbation. (C) Enriched cis-regulatory motifs in cells with different perturbations. I-cisTarget was used to predict regulatory motifs in the significantly (Padj <0.05) differentially expressed genes in cells with different gRNA, with indication of the normalized enrichment scores (NES). (D) Cell-cycle phase distribution per perturbation. Based on the normalized expression levels of cell-cycle marker genes, cells were classified in S, G2M or G1 phase and perturbations were ordered based on the fraction of cells in G1 phase. (E) Enrichment analyses showing enrichment scores of pro-T signature gene sets in the ranked list of differentially expressed genes in cells with Stat5b, Notch1 and Myc perturbation versus non-targeting control cells. The NES and false discovery rate q-value are indicated. gRNA: guide RNA; FDR: false discovery rate; NT: non-targeting control.

Figure 3.Combined analysis of single-cell and bulk transcriptome of cells with Spi1 perturbation. (A) Cell densities of BFP-positive cells (mean ± standard deviation [SD], N=3 replicates) over time for pro-T cells with Spi1 guide RNA (gRNA), compared with two non-targeting (NT) control gRNA. Data were analyzed by one-way analysis of variance at day 22. The statistical significance of P values is indicated by asterisks (***P<0.001, **P<0.01), or NS (not significant). (Spi1_gRNA1 vs. Spi1_gRNA3: P=0.9989; Spi1_gR-NA1 vs. NT1_gRNA1: P=0.0026; Spi1_gRNA1 vs. NT1_gRNA5: P=0.0010; Spi1_gRNA3 vs. NT1_gRNA1: P=0.0022; Spi1_gRNA3 vs. NT1_gR-NA5: P=0.0009; NT1_gRNA1 vs. NT1_gRNA5: P=0.8374). (B) Median fluorescence intensities (MFI) of SPI1 protein levels in cells with Spi1 gRNA versus NT by intracellular flow cytometry. Normalized MFI were calculated for AlexFluor568 emission. Three biological replicates were performed, with mean and SD shown. Comparison between groups was performed by a two-sided unpaired Student t test, P<0.0001. (C) Cell densities of GFP-positive cells (mean ± SD, N=3 replicates) over time for pro-T cells with SPI1 overexpression and an empty vector. Data were compared at day 15 by a two-sided unpaired Student t test, P=0.9662. (D) Overlap between significantly (Padj <0.05) differentially expressed genes in the bulk and single-cell data for Spi1 perturbed cells. P values for the overlaps were calculated using the hypergeometric probability distribution. (E) Gene set enrichment analysis showing normalized enrichment scores for the Hallmark gene sets in the ranked list of differentially expressed genes in cells with Spi1 gRNA versus NT control cells. (F) Enrichment analysis for the MYC gene expression signature in the single-cell dataset (upper panel) and the bulk dataset (lower panel) of Spi1 perturbation. The NES and false discovery rate q-value are indicated on the graphs. (G) Heat maps showing expression levels of genes within SCENIC regulons, in cells with Spi1 gRNA versus NT control cells. Small gene regulatory networks were determined based on co-expression patterns and transcription factor motifs. BFP: blue fluorescent protein; GFP: green fluorescent protein; OE: overexpression; EV: empty vector; RNA-seq: RNA sequencing; FDR: false discovery rate.

Characterizing the transcriptomic impact of recurrent mutations in T-cell acute lymphoblastic leukemia

We next selected 42 poorly characterized genes that are recurrently mutated in T-ALL patients with predicted loss-of-function alterations.^2,3,32,33 The CRISPR library was designed with five gRNA per gene, 11 NT gRNA and ten positive controls, resulting in a total of 231 gRNA (Online Supplementary Table S3). Single cells were sequenced at three timepoints (3, 7 and 14 days after transduction), generating a dataset of 18,032, 19,718 and 16,135 cells passing filtering at each respective timepoint, with a single gRNA detected in 5,280, 8,123 and 6,084 cells, respectively (Online Supplementary Figure S4A). Other quality metrics were comparable to those of the first CRISPR screen (Online Supplementary Figure S4B-D). For the majority of perturbations, downregulation of the target gene was observed, indicating efficient perturbation (Online Supplementary Figures S2 and S4E). In general, we aimed to investigate whether the transcriptomic impact of different perturbations converges on common downstream pathways. The data of both CRISPR screens were merged, with regression of batch and cell-cycle effects. Enrichment analyses were performed based on comparison of the pseudo-bulk transcriptomes using the previously defined signatures as reference (i.e., NOTCH, STAT, E2F and MYC)⁷ (Figure 4A, B, Online Supplementary Figure S4F). We observed a large group of perturbations (cluster A) with upregulation of both STAT and NOTCH, and downregulation of the E2F signature. Within this cluster, the effect on MYC transcriptional activity was variable, but often resulted in downregulation of the MYC signature. Similar to direct perturbation of Myc (Online Supplementary Figure S4F), these perturbations were commonly associated with upregulation of NOTCH and STAT signatures. Cluster A contained multiple genes for which their function in T-ALL remains unclear, such as Tspyl2, Trrap, Zbtb7a and Pbrm1 (Figure 4A). These genes clustered alongside multiple known epigenetic regulators, including Ctcf, Arid1a, Smarca4 and Ezh2, suggesting that these factors may impact the transcriptome of pro-T cells in a similar way.

By quantifying the number of cells over time, we analyzed the impact of each gRNA on cell proliferation (Figure 4C, Online Supplementary Figure S4G). While many of the perturbations had a negative effect on cell growth, only a few had a positive effect. Good correspondence was found between cell proliferation and cell-cycle phase as predicted by the single-cell RNA-sequening data (Figure 4D). Loss of Dnm2, Med12, Gata3, Smarca4, Arid1a, Ncor1 and Ctcf caused the strongest growth disadvantage, identifying these genes as dependency factors in pro-T cells. Data from the DepMap cancer dependency portal³⁴ also showed negative CRISPR scores for those genes in lymphoid cell lines (Figure 4E, Online Supplementary Figure S5). Using our single-cell data, we could link inactivation of these genes with a strong negative effect on MYC and E2F signatures (Figure 4A). In contrast, we observed a growth advantage for cells with Pten, Cdkn2a, Runx1, Bcl11b, Crebbp, Bcor, Setd2, Dnmt3a and Trp53 gRNA (Figure 4C). These genes were thus identified as tumor suppressors in pro-T cells and many are indeed known as such.^25,35-37 Inactivation of these tumor suppressor genes was associated with upregulation of either the MYC/E2F signatures (e.g., Pten, Runx1, Cdkn2a, Ptprc, Jun) or the STAT signature (e.g., Crebbp, Trp53, Bcl11b, Setd2, Dnmt3a) (Figure 4A, B).

We next applied SCENIC to study regulons associated with the different perturbations. We observed considerable variation in the expression of regulons of several ETS family members, such as Spi1, Fli1, Ets1 and Erg. Most strikingly, the Ets1 regulon, as well as the Ets1 gene, were upregulated specifically in cells with Spi1 gRNA and not in cells with any of the other perturbations (Figure 4F, Online Supplementary Figure S4H and S6). The Spi1 regulon was downregulated in cells with Spi1 and Pten gRNA, and upregulated in cells with Ptpn2, Gata3, Crebbp and Arid1a loss. We observed that the Cebpb regulon was active in about half of the perturbations, while the Cebpg regulon was active in roughly the other half of the cells, suggesting an antagonistic function (Online Supplementary Figure S6). Interestingly, we found that the Erg regulon was up-regulated in cells with a growth disadvantage (i.e., Arid1a, Ctcf, Dnm2, Kit, Kmt2c, Kmt2d, Myc and Spi1 inactivation) (Figure 4F). Upregulation of Erg seems to interfere with cell proliferation, and it was previously described that Erg knockout leads to a highly proliferative phenotype.¹⁷

Bcl11b is a tumor suppressor and cooperates with activated JAK/STAT signaling

We next focused on Bcl11b, which was identified as one of the strongest tumor suppressor genes in pro-T cells in our CRISPR screen (Figure 4C). To validate the results of the screen, we performed independent Bcl11b inactivation by two gRNA, confirming the increase in cell proliferation (Figure 5A) and the direct effect on BCL11B protein level (Figure 5B). Similarly, Bcl11b knockdown by shRNA also resulted in increased cell proliferation, although at a more modest rate compared to that caused by CRISPR perturbation (Online Supplementary Figure S7A). In contrast, Bcl11b overexpression significantly impaired cell growth, confirming the growth-suppressing role of BCL11B (Figure 5C).

Figure 4.Perturbations can be clustered based on their gene expression profiles. (A) Heatmap of pro-T gene signature clusters NOTCH, STAT, E2F and MYC. The data of both CRISPR screens were merged and the normalized enrichment scores (NES) were calculated for each perturbation. Unsupervised hierarchical clustering was used to cluster perturbations based on similarity of the NES. For the blank boxes, no NES could be calculated. (B) Volcano plots of the gene sets E2F, MYC, NOTCH and STAT for cells with perturbation of Bcl11b, Ctcf, Spi1 and Arid1a. Dotted lines are thresholds for significantly (Padj <0.05) differentially expressed (|log2Fc|>0.5) genes. (C) Evolution of the fraction of cells per perturbation over time at days 3, 7 and 14 of screen 2. The Y-axis shows data on a logarithmic scale. Non-targeting controls (NT1, NT2, and NT3) controls are indicated in gray font. Perturbations are ordered based on the number of cells at day 14. (D) Cell-cycle phase distribution per perturbation in screen 2. Based on normalized expression levels of cell-cycle marker genes, cells were classified in S, G2M or G1 phase and perturbations were ordered based on the fraction of cells in G1 phase. (E) Average CRISPR scores from the DepMap database for all genes in lymphoid leukemia and lymphoma cell lines. Negative scores mean a negative impact on cell proliferation, while genes with positive scores cause a proliferative advantage. Genes on the X-axis were ordered by increasing gene effect score. (F) Enrichment of SCENIC regulons per perturbation. Stacked violin plot displaying the distribution of AUCell values for each regulon and each perturbation. AUCell values represent enrichment of each predicted regulon in each single cell. The color intensity is a measure of the median expression of each regulon in each perturbation. This figure shows a selection of perturbations and regulons, with complete data displayed in Online Supplementary Figure S6.

Single-cell RNA sequencing indicated that loss of Bcl11b in pro-T cells was associated with increased STAT and E2F signatures (Figure 4A, B). For in-depth characterization of the impact of Bcl11b inactivation, bulk RNA sequencing was performed. We found 2,690 significantly differentially expressed genes, with significant overlap with the single-cell dataset (Figure 5D, Online Supplementary Figure S7B, Online Supplementary Table S2). These data were integrated with a publicly available chromatin immunoprecipitation-sequencing dataset to define high-confidence direct BCL11B target genes (Figure 5E).³⁸ Enrichment analysis of these direct target genes demonstrated significant upregulation of inflammatory signatures, NF-kB and STAT signaling (Figure 5F). Cis-regulatory motif analysis showed that RELA, STAT2 and SPI motifs were positively enriched, while EZH2, FOXO1 and SOX11 motifs were negatively enriched (Figure 5G). In addition to enrichment of the NF-kB pathway and RELA motif, we found increased P65 protein level and phosphorylation, as well as elevated expression of the Rela and Nfkb2 regulons (Figure 5G, H, Online Supplementary Figure S7C). Additionally, several genes of the BCL2 family (Bcl2a1, Bim, Mcl1) were upregulated (Online Supplementary Table S2), which are negative regulators of apoptosis and known NF-kB target genes. Besides NF-kB activation, the data showed significant upregulation of STAT signaling induced by Bcl11b loss (Figure 5F-H). Together, these data identified Bcl11b as a strong tumor suppressor in pro-T cells, associated with deregulated NF-kB and STAT signaling.

To validate these findings in samples from patients with T-ALL, we used a publicly available RNA-sequencing dataset of 264 T-ALL samples.² We divided the patients into three groups based on BCL11B expression levels (Figure 6A, Online Supplementary Figure S8) and found that T-ALL cases with the lowest BCL11B expression also showed significant enrichment of NF-kB and STAT signaling (Figure 6B). In concordance, we found significantly reduced BCL11B expression in T-ALL patients with alterations in the JAK/STAT signaling pathway (Figure 6C). This two-way correlation led us to hypothesize that BCL11B inactivation could have a synergistic effect in combination with activating JAK/STAT alterations. Indeed, Bcl11b perturbation in combination with overexpression of the activating JAK3^M511I mutant was able to drive cytokine-independent growth of pro-T cells, while single mutants did not have transforming potential (Figure 6D). Treatment with ruxolitinib demonstrated that pro-T cells with loss of Bcl11b had decreased sensitivity towards pharmacological JAK targeting (Figure 6E), in agreement with a direct effect of BCL11B on STAT transcriptional activity.

Discussion

This work describes a comprehensive study of a number of recurrently mutated genes in T-ALL, which typically harbor loss-of-function mutations but remain poorly characterized. The aim was to determine whether these genes have tumor suppressive characteristics and to correlate this with specific effects on gene expression. Single-cell CRISPR screening enabled us to study the effects of gene perturbation on proliferation as well as transcriptional regulation. An ex vivo model of primary pro-T cells was used, as these immature cells are likely closely related to the cells of origin of T-ALL. We initially validated the methodology by targeting well-characterized T-cell regulators, achieving efficient downregulation of the majority of genes, as well as significant differential expression of known downstream target genes, pathways and enriched cis-regulatory motifs. Furthermore, these data enabled us to define the core transcriptional signatures driving pro-T-cell growth.

We subsequently performed a CRISPR screen targeting 42 genes that are frequently altered in T-ALL patients. By assessing cell number over time, these perturbations were linked with their effects on cell proliferation. Interestingly, we observed upregulation of the STAT signature by the majority of perturbations. It is well-established that deregulated JAK/STAT signaling is an important driver of T-ALL, as illustrated by the numerous mutations in this pathway (in JAK3, JAK1, IL7R, PTPN2, etc.) that are detected in T-ALL patients.² Based on our CRISPR data, we could hypothesize that many genetic mutations eventually lead to deregulated STAT signaling, thereby contributing to disease development. Moreover, we identified a large cluster of genes exhibiting upregulation of both the STAT and NOTCH signatures, with a smaller subcluster characterized by additional downregulation of MYC and associated with reduced cell proliferation. As MYC is a master transcriptional regulator in T cells, it is not surprising that its downregulation is associated with strong fitness effects. For instance, our data showed a strong growth disadvantage for cells with Ctcf inactivation, coinciding with MYC downregulation. This is in line with previous data illustrating a role for CTCF in chromatin looping of the MYC enhancer to the MYC promoter.³⁹ Remarkably, many of the genes in cluster A encode epigenetic regulators, including Arid1a and Smarca4 as part of the SWI/SNF complex and Kdm6a, Kmt2c and Asxl2 as regulators of histone (de)methylation.^40-42 Our data suggest that other genes in this cluster, such as Tspyl2 or Zbtb7a, may have similar transcriptional effects. In addition to genes in cluster A, we identified Spi1 as an essential gene for cell survival and linked this dependency to MYC. A relationship between SPI1 and MYC was previously described in myeloid cells in which Spi1 perturbation leads to downregulation of MYC targets due to loss of SPI1 binding to MYC enhancers.⁴³ Additionally, we observed upregulation of Ets1, as well as regulons associated with other ETS factors. These data complement work from Rothenberg’s laboratory who previously characterized the role of several transcription factors, including PU.1, in mouse T-cell development.¹⁷

Figure 5.Bcl11b inactivation drives cell proliferation and upregulation of NF-kB and STAT signaling. (A) Cell densities of BFP-positive cells (mean ± standard deviation [SD], N=3 replicates) over time for pro-T cells with Bcl11b guide RNA (gRNA), compared with two non-targeting (NT) control gRNA. Data were analyzed by one-way analysis of variance at day 22. The statistical significance of P values is indicated by asterisks (**P<0.01, *P<0.05), or NS (not significant). (NT1_gRNA1 vs. NT1_gRNA5: P=0.9993; NT1_gRNA1 vs. Bcl11b_gRNA1: P=0.0037; NT1_gRNA1 vs. Bcl11b_gRNA5: P=0.0367; NT1_gRNA5 vs. Bcl11b_gRNA1: P=0.0043; NT1_gRNA5 vs. Bcl11b_gRNA5: P=0.0434; Bcl11b_gRNA1 vs. Bcl11b_gRNA5: P=0.3575). (B) Protein levels of BCL11B as measured by intracellular flow cytometry. Normalized median fluorescence intensities were calculated for AlexaFluor555 emission. Three biological replicates were performed, with mean and SD shown. Comparisons between groups were performed by a two-sided unpaired Student t test, ****P<0.0001. (C) Cell densities of GFP-positive cells (mean ± SD, N=3 replicates) over time for pro-T cells with Bcl11b over-expression and an empty vector. Data were compared at day 10 by a two-sided unpaired Student t test, P=0.0004. (D) Overlap between significantly (Padj <0.05) differentially expressed genes in the bulk and single-cell data for cells with Bcl11b inactivation compared to NT cells. P values for the overlaps were calculated using the hypergeometric probability distribution. (E) Overlap between chromatin immunoprecipitation-sequencing peaks obtained from a publicly available dataset of BCL11B binding regions, and the significantly differentially expressed genes in the bulk RNA-sequencing data. P values for the overlaps were calculated using the hypergeometric probability distribution. (F) Gene set enrichment analysis on the set of direct Bcl11b target genes. Normalized enrichment scores and P values are shown for the Hallmark gene sets. (G) Predicted cis-regulatory binding motifs in the significantly (Padj <0.05) up- and downregulated genes in cells with Bcl11b inactivation. (H) Heatmaps showing expression levels of genes within SCENIC regulons in cells with Bcl11b gRNA versus NT control cells. BFP: blue fluorescent protein; MFI: median fluorescence intensity; GFP: green fluorescent protein; EV: empty vector; OE: overexpression; RNA-seq: RNA sequencing; ChIP-seq: chromatin immunoprecipitation sequencing; NES: normalized enrichment score.

Figure 6.Bcl11b levels are correlated with STAT pathway activation and cause transformation of pro-T cells. (A) Normalized expression level of BCL11B in publicly available RNA-sequencing data of patients with T-cell acute lymphoblastic leukemia (T-ALL). Patients were subdivided into three groups based on expression level of BCL11B (i.e., average – high – low). (B) Gene set enrichment analysis comparing patients with low BCL11B expression with those with high expression level. Normalized enrichment scores and P values are shown for the Hallmark gene sets. (C) Normalized expression level of BCL11B in publicly available RNA-sequencing data of T-ALL patients with alterations in the JAK/STAT signaling pathway compared to patients without mutations in this pathway. The P value was calculated with an unpaired Student t test. (D) Cell densities (mean ± standard deviation, N=3 replicates) of pro-T cells over time with the following conditions: JAK3^M511I + non-targeting control (NT), JAK3^M511I + gBcl11b, empty vector (EV) + gBcl11b and EV + NT. The left panel shows cells cultured in the presence of IL7 and the right panel shows IL7-deprived cells. Data were analyzed by one-way analysis of variance at day 14 (left panel) or day 7 (right panel). The statistical significance of P values is indicated by asterisks (****P<0.0001, ***P<0.001, **P<0.01, *P<0.05) or NS (not significant). Left panel: JAK3^M511I + NT vs. JAK3^M511I + Bcl11b gRNA: P=0.0063; JAK3^M511I + NT vs. EV + NT: P=0.9981; JAK3^M511I + NT vs. EV + Bcl11b gRNA: P=0.0035; JAK3^M511I + Bcl11b gRNA vs. EV + NT: P=0.0078; JAK3^M511I + Bcl11b gRNA vs. EV + Bcl11b gRNA: P=0.9634; EV + NT vs. EV + Bcl11b gRNA: P=0.0043. Right panel: JAK3^M511I + NT vs. JAK3^M511I + Bcl11b gRNA: P<0.0001; JAK3^M511I + Bcl11b gRNA vs. EV + NT: P<0.0001; JAK3^M511I + Bcl11b gRNA vs. EV + Bcl11b gRNA: P<0.0001. (E) Dose-response curve for pro-T cells with JAK3^M511I and Bcl11b gRNA treated with ruxolitinib phosphate for 60 h. Relative viability was calculated by normalizing the data based on the dimethylsulfoxide condition of each group. NES: normalized enrichment score; mut: mutated; WT: wild-type; IL7: interleukin-7.

Several of the genes in our screen were here identified as essential genes in pro-T cells, despite harboring loss-of-function mutations in T-ALL patients. Our findings are in line with CRISPR data from the DepMap database in which many of those genes were similarly recognized as essential in lymphoid leukemia and lymphoma cell lines.³⁴ Many of these genes thus seem to behave as haploinsufficient tumor suppressor genes in T-ALL, with complete inactivation by CRISPR not being tolerated. This is likely explained by the fact that many loss-of-function mutations in T-ALL patients are either heterozygous or are point mutations leading to specific changes in protein function, whereas CRISPR typically causes complete knock-out. Still, we demonstrated that mild knockdown by shRNA also had similar effects as CRISPR perturbation on pro-T-cell proliferation. It was previously shown that T-ALL patients have a complex genetic landscape and can carry up to 20 different mutations, with many of these alterations showing either significant co-occurrence or mutual exclusivity.^2,44 Single mutations are often not sufficient to drive leukemia development in mouse models and multiple cooperating mutations are required to initiate disease.^6-8,10,45,46 We hypothesize that this genetic complexity is an important factor that shapes the disease. As it was previously shown that the order of mutational acquisition in ALL is not random,⁴⁷, ^{4 8} we reason that the timing of each mutation is crucial and some phenotypes may only manifest on a background in which earlier driving mutations are present.

We identified Bcl11b as a strong tumor suppressor gene in pro-T cells and demonstrated a direct connection between effects on cell proliferation and gene expression. There are conflicting data on the role of BCL11B which most likely reflects the results obtained in various model systems using various technologies. Our CRISPR screen allowed us to study the inactivation of Bcl11b in primary immature T cells. In physiological conditions, Bcl11b expression initiates around the commitment phase and pushes progenitor cells forward through T-cell development.⁴⁹ BCL11B insufficiency has been described to lead to loss of commitment, reduced sensitivity to NOTCH stimulation and allows reprogramming of T cells into natural killer or myeloid fates.^49,50 In T-ALL, BCL11B is deleted in 3% of patients and shows mutations in around 10% of cases, often co-occurring with aberrant TLX1 expression. In other cases, a t(5;14)(q35;q32) translocation juxtaposes the BCL11B regulatory region to the TLX3 coding sequence, thereby inactivating BCL11B.^49,51 In our dataset, as well as publicly available RNA-sequencing data from T-ALL patients, significant upregulation of the NF-kB and STAT pathways was observed in cells with low Bcl11b levels. Moreover, we also found low expression of BCL11B in T-ALL patients’ samples with JAK/STAT pathway alterations, indicating that reduced levels of BCL11B might cooperate with JAK/STAT pathway activation. Activating JAK3 mutations are frequent events in T-ALL and it was previously demonstrated that JAK3 signaling can cooperate with other factors, such as Suz12 inactivation and Il7r mutations.^3,6,9,46,52 In the current study, we found that inactivation of Bcl11b enhanced the oncogenic effects of mutant JAK3 and drove transformation to cytokine-independent growth. In addition, cells with Bcl11b inactivation showed decreased sensitivity to JAK inhibition, irrespective of JAK3 mutant status. These data, in combination with the single-cell RNA-sequencing data, indicate that BCL11B functions as an important negative regulator of STAT transcriptional activity in immature T cells. The absence of Bcl11b results in increased activation of STAT transcriptional activity, making the cells more vulnerable to JAK pathway mutations, yet simultaneously reducing their sensitivity to JAK kinase inhibitors.

Our data demonstrate that single-cell CRISPR screening is a powerful technique for parallel characterization of the transcriptional alterations resulting from numerous gene perturbations. In this way, we characterized Bcl11b as a potent tumor suppressor gene in T-ALL with a link to JAK/STAT and Spi1 as an essential gene associated with a direct effect on regulation of Myc. These data offer insight into the function of uncharacterized genes and can be linked with other characteristics such as effects on proliferation or cell cycle.

Footnotes

• Received December 27, 2023
• Accepted May 17, 2024

Correspondence

J. Cools jan.cools@kuleuven.be

S. Demeyer sofie.demeyer@kuleuven.be

Disclosures

No conflicts of interest to disclose.

Funding

This work was funded by FWO (11L9122N) and KU Leu-ven PhD fellowships to SM, a Foundation Against Cancer postdoctoral fellowship to SD, VIB Tech Watch Technology Implementation fund (to JC) and KU Leuven funding (C14/18/104, to JC)

Acknowledgments

This study was aided by the KU Leuven Genomics Core facility and VIB Flow Core facility.