Single-cell analysis identifies dynamic gene expression networks that govern B cell development and transformation

Abstract

Integration of external signals and B-lymphoid transcription factor activities organise B cell lineage commitment through alternating cycles of proliferation and differentiation, producing a diverse repertoire of mature B cells. We use single-cell transcriptomics/proteomics to identify differentially expressed gene networks across B cell development and correlate these networks with subtypes of B cell leukemia. Here we show unique transcriptional signatures that refine the pre-B cell expansion stages into pre-BCR-dependent and pre-BCR-independent proliferative phases. These changes correlate with reciprocal changes in expression of the transcription factor EBF1 and the RNA binding protein YBX3, that are defining features of the pre-BCR-dependent stage. Using pseudotime analysis, we further characterize the expression kinetics of different biological modalities across B cell development, including transcription factors, cytokines, chemokines, and their associated receptors. Our findings demonstrate the underlying heterogeneity of developing B cells and characterise developmental nodes linked to B cell transformation.

Fig. 1. Identifying B-cell development stages using scRNAseq and CITE-Seq.

a Schematic of the experimental setup. WT wild type. b UMAP dimension-reduction projection of all cells (n = 7454) from two wild-type C57BL6 mice. Fourteen clusters were identified and the corresponding population names for each cluster are listed (Pre-BCRd pre-BCR-dependent, Pre-BCRi pre-BCR-independent) c Feature plot of cells that are labeled according to their cell cycle status based on gene expression (left) and Mki67 transcription expression (right). Cell cycle status was determined by the average expression of gene sets representing each cell cycle, including postmitotic G1 phase (G1PM), S phase, or G2/M phase, and cells that did not harbor any of these signatures were labeled as G0/G1. Color scale represents natural log-transformed SCTransform corrected counts. d Feature plot of cells for their CITE-Seq/ADT antibody expression (top row) and the corresponding gene transcript expression (bottom row). Color scale represents centered natural log transformation across cells. adt antibody-derived tags.

Fig. 2. Transcriptional signatures of pre–pro B cells and committed pro-B cells.

a Highlighted populations for pre–pro B and cycling pro-B cells (left). Feature plots for early B-cell markers (right). b Heatmap of genes that are differentially expressed between the pre–pro B cells and the cycling pro-B cells. Scale represents normalized counts centered and scaled across cells. c Feature plot of EBF1-target genes. d Differentially upregulated genes in the CD93 ADT high versus CD93 ADT low pre–pro B-cell cluster. CD93 high versus low was based on a median split. No genes were significantly downregulated in the CD93 high population. e Feature plot of genes denoting the pro-B-cell stage (Kit), or genes uniquely expressed during pro-B cells. f Violin plot of the expression of surrogate light chains, Vpreb1 and Igll1 across the B-cell development stages. Color scale in a, c, e, and f represents natural log-transformed SCTransform-corrected counts.

Fig. 3. The transcriptome of the pre-BCR-dependent expansion stage is governed by MYC-associated gene expression networks and requires repression of Ebf1 expression.

a Highlighted population of the pre-BCR-dependent cluster (left) and feature plots for the pre-BCR-dependent markers highly expressed during this stage (right). b Violin plot of B-lymphoid transcription-factor expression, including Ebf1, Pax5, and Ikzf1 across B-cell development. c Violin plot of additional genes that are downregulated during pre-BCR-dependent proliferation, which includes Il7r, and EBF1-target genes such as Cd79a and Cd79b. d Landscape In silico Deletion Analysis (LISA) to predict the transcriptional regulator of the top 100 differentially upregulated genes during pre-BCR-dependent proliferation (left table). Gene-set enrichment analysis was performed to identify gene sets from the Molecular Signature Database that were enriched in the pre-BCR-dependent cells. Comparison between the pre-BCR-independent cells (left side) and the pre-BCR-dependent cells was performed (right side). Statistical significance was calculated using one-sided Wilcoxon rank-sum test. e EBF1-binding sites at the promoters of Nrgn, Slc7a5 and Slc3a2 with indicated MACS peak calls. Data were obtained from GSM2863146. f Log2-transformed FPKM expression values obtained from RNA-seq of wild-type and Pax5^+/− × Ebf1^+/− (PE) leukemic progenitor B cells. Progenitor B cells were obtained via negative selection of CD11c, TER119, Ly6G, Ig Kappa, and Ig Lambda and positive selection of CD19. Statistical significance was determined using a two-tailed unpaired student t-test for Ybx3 (P = 0.0024) and Slc7a5 (P = 0.0001). RNA-seq data were obtained from GSE148680. A two-tailed Mann–Whitney test was performed for Slc3a2 (P = 0.0083) due to non-normal distribution. n = 7 biologically independent samples over one independent experiment. Color scales in a, b, and c represent natural log-transformed SCTransform-corrected counts.

Fig. 4. Pre-BCR-independent proliferation is distinct from pre-BCR-dependent proliferation.

a Highlighted pre-BCRi I and pre-BCRi II populations (left). Violin plots of expression of ADT-CD43 and ADT-CD25 across the B-cell development stages (right). b Feature plot and violin plot for Cd74, Cd44, and Bach2. c Heatmap of differentially expressed genes between the pre-BCR-dependent (Pre-BCRd), pre-BCR-independent I (Pre-BCRi I), and pre-BCR-independent II (Pre-BCRi II) S, or G2/M phase. Scale represents normalized counts centered and scaled across cells. d Volcano plot showing differentially regulated genes between the pre-BCRd cluster and the pre-BCRi I cluster. e Differentially regulated genes between the pre-BCRi I cluster and the pre-BCRi II cluster. Color scales in a, b represent natural log-transformed SCTransform-corrected counts.

Fig. 5. B-cell differentiation and maturation.

a Highlighted B-cell clusters undergoing V(D)J recombination (left). Feature plot of genes involved in V(D)J recombination or specific for pre-B-cell expression (right). b Highlighted B-cell clusters for late B-cell maturation (left). Feature plot of genes highly expressed in immature or mature B cells (right). Color scale in a,b represents normalized counts centered and scaled across cells.

Fig. 6. Pseudotime analysis illustrates the kinetics of transcription-factor expression and gene modules that are differentially expressed across the B-cell development trajectory.

a Cell cycle-related genes were regressed within Monocle3 and UMAP dimension reduction was performed. Clustering was performed based on the Monocle3 clusters (left). Cells from Monocle-based clustering labeled with Seurat cluster labels (right). b Pseudotime values were calculated and plotted. c Module analysis to demonstrate gene modules that change across the B-cell developmental trajectory. A total of 33 modules and their expression intensity for each stage are shown. Color scale represents normalized module gene expression. d Gene ontology term analysis of selected modules 1, 6, and 8. e Expression of B-lymphoid transcription factors and epigenetic factors across B-cell development stages. f Expression of cytokine, chemokine and cytokine/chemokine receptors across the B-cell developmental trajectory.

Fig. 7. B-cell developmental gene networks correlate with various B-ALL subtypes and associate with prognosis.

a Heatmap of cluster gene sets and various human B-ALL subtypes. n_number denotes the number of patient samples falling in each B-ALL subtype. Average z-score of cluster marker genes was calculated for each cluster and hierarchical clustering was performed for both B-ALL subtype and average cluster z-score. b Monte-Carlo-based statistical testing of the hierarchical clustering was performed for row (B-cell development stages) and column (human B-ALL subtypes). Significant p-values (<0.05) are shown and colored as a red line. c Overlapping differentially upregulated gene markers in the pre-BCR-dependent clusters from Seurat and Monocle are plotted (y axis). The expression of these genes was queried across various B-ALL subtypes (x axis). The number of B-ALL samples (n) is listed after each subtype. d Survival curve of pediatric B-ALL (left, time in years) and adult B-ALL (right, time in months). Low and high expression correspond to below-median and above-median expression of YBX3, respectively. Survival data were obtained from Prediction of Clinical Outcomes from Genomic Profiles (PRECOG) e Proposed model of B-cell development.