Count values generated from RNA-seq data of leukemia patients were downloaded from St. Jude Cloud (23). FPKM values were calculated from count values with an in-house perl script. Cases with GATA2 FPKM values >3-fold standard deviation of the mean in each subtype were categorized as GATA2 activation. Protein-coding genes with FPKM >1 in at least 30% of patients were defined as expressed genes and collected for further analysis. Batch effects introduced by different sequencing protocols (Supplementary Table 1) were corrected with DESeq2. Differentially expressed gene (DEG) analysis was performed on the expressed genes using DESeq2 (24) with subtype as a covariate. Genes with fold change >1.5 and p-value < 0.05 were identified as DEGs.

Gene set enrichment analysis of genes up-regulated and down-regulated upon GATA2 activation was performed using ToppGene (25). The categories with top significance enriched in Gene Ontology (GO) terms and ToppCell Altas which points to cell type in ToppGene gene set were shown in figures.

GATA2 ChIP-seq data (GSM1600544) from the LNCaP cell line was used to calculate the distance between GATA2 peaks and the transcription start site (TSS) of the nearest gene to determine the genomic range for motif analysis. According to the results, the 200 bp sequences centered on the TSS of protein coding genes across the whole genome were analyzed with FIMO (26) for the presence of the transcription factor binding motif of the GATA family. Predictions with a p-value < 0.0001 were considered positive.

The heterozygous markers and RNA-seq BAM file for SJALL043839_D1 were obtained from previously published data (27). Allelic specific expression (ASE) analysis was performed with cis-X (28) following the user's instructions. The ASE of GATA2 was calculated with in-house R script.

The protein structure of the in-frame KMT2A-USP2 fusion was visualized with ProteinPaint (29). The topologically associating domain (TAD) structure from high-through chromosome conformation capture (Hi-C) data from the Nalm6 and Kelly cell line, and the coverage plot of whole-genome sequencing (WGS) data of tumor and matched remission sample of SJALL043839 were visualized with GenomePaint (30).

The DEG analysis was performed with DESeq2. GATA family member expression between different lineages of leukemia was compared using Wilcoxon test. The correlation in expression between GATA2 and potential target genes was calculated with Pearson correlation. All analysis were performed in R (version 3.5.3).

We examined the RNA expression of three hematopoietic GATA transcription factors (GATA1, GATA2, and GATA3) in a total of 1,644 pediatric leukemia patients with AML (n = 297), T-ALL (n = 99), and B-ALL (n = 1,248) using RNA-seq data available at St. Jude Cloud (23) (Figure 1A). Consistent with their established roles in hemopoiesis (16, 31), the three GATA transcription factors showed lineage specific patterns of transcription. GATA1 and GATA2 were more highly transcribed in AML as compared with B-ALL (p-value < 0.0001, Wilcoxon test) and T-ALL (p-value < 0.0001, Wilcoxon test), and GATA3 showed the highest level of transcription in T-ALL (p-values < 0.0001 in comparison with AML and B-ALL, Wilcoxon test). On the other hand, transcription levels of all three GATA genes were low in B-ALL. In particular, we found that GATA2 was not transcribed in the majority of B-ALL patients, with FPKM values <1 in 79.4% (991 out of 1,248) cases (median FPKM = 0.297). The lack of GATA2 transcription in B-ALL was in accordance with its role in common lymphoid progenitor cells during normal hemopoiesis (32). However, it was noticeable that GATA2 was actively transcribed in a subset of B-ALL patients (n = 257 with FPKM > 1).

We next asked if GATA2 transcription was enriched in a specific B-ALL subtype. Toward this end, we checked transcription levels of GATA2 in 629 of the 1,248 B-ALL cases mentioned above, for which definitive subtype information was available. First of all, we found that active transcription of GATA2 was associated with ETV6-RUNX1 rearranged B-ALL (ETV, Figure 1B, median FPKM = 2.26, p < 0.0001, Wilcoxon test, ETV compared with non-ETV). Secondly, aberrantly high GATA2 transcription was observed in a subset of B-ALL patients across different subtypes. We next identified the cases with outlier GATA2 transcription defined as GATA2 transcription levels greater than three standard deviations of the mean in each subtype. Aberrant transcription was analyzed within each subtype to minimize the variation across different subtypes. According to this criterion, a total of 13 B-ALL cases with aberrant activation of GATA2 were identified from 7 subtypes (Figure 1B), accounting for 2.07% of B-ALL cases. These 13 B-ALL cases included five with high hyperdiploid (HYPER), two with KMT2A rearrangements (MLL), two with ETV, and one each with DUX4 rearrangement (DUX4), low hypodiploid (HYPO), BCR-ABL1 (PH), and BCR-ABL1-like (PHL) (Figure 1C). In addition, no outlier GATA2 transcription was identified for the following four subtypes: TCF3-PBX1 (E2A), MEF2D rearrangements (MEF2D), PAX5 alterations or mutations (PAX5), and ZNF384 rearrangements (ZNF384) (Figure 1B; Supplementary Table 1). To further investigate the effect of GATA2 activation, we selected 584 patients from seven subtypes with at least one GATA2 outlier case identified. Among these 584 cases, 13 cases with aberrant activation of GATA2 were classified as the GATA2-outlier group and the remaining cases (n = 571) as the GATA2-normal group.

To understand the impact of dysregulated GATA2 transcription on B-ALL, we performed differential expression analysis between the GATA2-outlier group and the GATA2-normal group using DESeq2 with B-ALL subtypes as a covariate. Out of 10,505 expressed protein coding genes (FPKM >1 in >30% patients), we identified 1,150 differentially expressed genes (DEGs), including 699 up-regulated and 451 down-regulated, in GATA2-outlier cases (p-value < 0.05 and |fold change| >1.5; Figure 2A; Supplementary Table 2). We firstly checked the transcription of B cell and myeloid lineage specific markers commonly tested in clinical diagnosis in the DEG result (Supplementary Table 3). Seven out of twelve myeloid markers used in clinical test were significantly up regulated in GATA2-outlier cases (Fisher's exact test, p < 0.001) (Figures 2A,B), while none of the 12 myeloid markers were down-regulated. On the other hand, B-lineage markers were significantly down regulated (Fisher's exact test, p < 0.001) (Figure 2A; Supplementary Figure 1A) and none of these B-lineage markers showed up-regulated expression. We next performed gene set enrichment analyses on up-regulated and down-regulated DEGs separately using ToppGene tool (25). Consistently, the 699 up-regulated DEGs were significantly enriched in myeloid leukocyte-associated enrichment categories including both GO (terms with a p-value < 10⁻³³ are listed in Figure 2C) and ToppCell Atlas (categories with a p-value < 10⁻³³ are listed in Figure 2D), suggesting the activation of a myeloid lineage-related transcription circuit following GATA2 dysregulation (15). In contrast, the 451 down-regulated DEGs were significantly enriched for the early B-lineage precursor cells stage (Supplementary Figure 1B) from which B-ALL leukemic cells originate. These observations indicated a potential switch from B-cell lymphocytes related transcription circuit to myeloid lineage related transcription circuit upon GATA2 activation. Down-regulated DEGs were also significantly enriched for cell cycle process (Supplementary Figure 1C), implying that GATA2 could be involved in the proliferation of blasts in GATA2-outlier patients.

To investigate the targets of GATA2 activation that are potentially involved in B-ALL, we combined genome-wide computational prediction of GATA family transcription factor binding motifs and the differential expression analysis described above. We first analyzed GATA2 chromatin immunoprecipitation sequencing (ChIP-seq) profiles of LNCaP cancer cell line published previously (33) to estimate the distance between each GATA2 binding site and its potential target gene. Only protein coding genes were included in this analysis. We found that among the GATA2 peaks with a nearby potential target gene (within 2 kb of the peak), 37.3% were located within ±100 bp of the transcription start site (TSS) of the nearby gene (Supplementary Figure 2A). This suggested that GATA2 targets downstream genes by predominantly binding to the ±100 bp of TSS region. Therefore, we scanned the TSS ±100 bp sequences of all protein coding genes across the human genome to predict GATA2 binding sites. This resulted in the identification of 313 genes with one or more GATA family binding motifs that may be direct targets of GATA factors (p-value < 0.0001 with FIMO, Supplementary Table 4). The predictions were confirmed by gene set enrichment analysis; the 313 genes showed significant enrichment in experimentally derived GATA binding sites from the ChIP-seq analysis in ENCODE project (34, 35) (Supplementary Figure 2B), supporting the prediction of GATA2 targets.

We next integrated the DEGs from differential expression analysis with the 313 targets of GATA factors from computational prediction. We identified 46 potential targets dysregulated by aberrant activation of GATA2, including 13 down-regulated DEGs (Supplementary Figure 2C) and 33 up-regulated DEGs (Figure 3A; Supplementary Table 2). We focused on the up-regulated DEGs in this analysis. Notably, there were several known tumor-associated genes (highlighted in Figure 2A) among these 33 target genes, including BAG3, EPOR, and KLF1. BAG3, which was reported to be an anti-apoptotic gene involved in leukemic cell survival and response to therapy (36), was among the most significant DEGs between GATA2-outlier and GATA2-normal B-ALL patients (p-value = 7.9 × 10⁻¹⁴; log2FC = 3.4, Figure 3B). And there was a significant positive correlation between GATA2 and BAG3 transcription in B-ALL (r = 0.35, p-value = 3.09e-18, Figure 3C). Furthermore, GATA2 binding was observed in ChIP-seq experiments from multiple tumor cell lines, namely LNCaP, HUVEC, SHSY5Y, and K562 (Figure 3D). Taken together, these observations support the hypothesis that BAG3 is a target of aberrantly activated GATA2 in B-ALL. Similar observations were made for EPOR and KLF1; both genes showed significant co-expression with GATA2 (Figure 3C) and positive GATA2 binding observed for multiple cancer cell lines in ChIP-seq experiments (Figure 3D). EPOR is a well-established B-ALL-associated gene whose abnormal activation in pre-B cells increases cell survival by activating the JAK-STAT pathway (37), whereas KLF1 is associated with poor survival in AML (38). The up-regulation of these genes by GATA2 activation may play a tumor-promoting role in leukemia and reinforce the potential function of GATA2 in leukemogenesis in B-ALL.

In the 13 GATA2-outlier B-ALL patients, two were MLL subtype (SJBALL021549_D1 and SJINF066_D). Notably, we found that both patients carried the same KMT2A-USP2 fusion. Further checking confirmed that these patients were the only two carrying this fusion in the cohort analyzed (Supplementary Table 1), indicating that GATA2 could be trans-activated by KMT2A-USP2. KMT2A-USP2 is a relatively rare type of KMT2A rearrangement and was recently reported as a recurrent fusion in pediatric B-ALL (39). In the two cases identified in this analysis, the breakpoint of the fusion was located within intron 22 of KMT2A (NM_005933) and intron 2 of USP2 (NM_004205), which resulted in an in-frame fusion consisting of exons 1–21 of KMT2A and exons 3–13 of USP2 (Supplementary Figure 3). The fusion protein contains the entire deubiquitination domain of USP2, indicating that USP2 might be the functional executor of this chimeric protein. Interestingly, wild-type USP2 was not transcribed in majority of B-ALL patients (Figure 4A) and could be activated through fusion with KMT2A (Figure 4B), with active transcription of USP2 occurring through hijacking of the KMT2A promoter (Supplementary Figure 3). In addition, we found high levels of wild-type USP2 transcription in couple of B-ALL patients with other MLL subtype (Figures 4A,B), indicating USP2 could also be aberrantly activated in KMT2A rearranged B-ALL. Notably, GATA2 outlier expression was found in patients with the KMT2A-USP2 fusion but not in other KMT2A rearranged patients with high expression of wild-type USP2 (Figure 4C), suggesting that GATA2 transcription was activated by the chimeric fusion protein as a whole.

Meanwhile, we were able to identify another MLL B-ALL case, SJALL043839_D1, with aberrant GATA2 transcription from a previously published ALL cohort (27), in which both RNA-seq and whole-genome sequencing (WGS) data were available for tumor cells. By integrating WGS and RNA-seq data with cis-X (28), we found significant allele-specific expression (ASE, Figure 5A) as well as outlier expression of GATA2 in this case compared with the remaining B-ALL patients from the same cohort (Figure 5B). Further analysis identified a somatically acquired 101 kb focal deletion located 285 kb downstream of the TSS of GATA2 in the tumor cells (Figure 5C). Notably, GATA2 was the only cis-activated gene in the genomic neighborhood surrounding this deletion (Supplementary Figures 4A–C), indicating that GATA2 was the target. Furthermore, in-situ high-through chromosome conformation capture (Hi-C) data for both the B-ALL cell line Nalm6 and neuroblastoma cell line Kelly (Figure 5C) from GenomePaint (30) showed that the deletion overlapped with the boundary of topologically associating domain (TAD) harboring GATA2. Disruption of the TAD structure by the deletion might aberrantly connect the GATA2 promoter to the active regulatory DNA elements in the adjacent TAD, thus contributing to GATA2 dysregulation in cis.