Differential Transcriptional Reprogramming by Wild Type and Lymphoma-Associated Mutant MYC Proteins as B-Cells Convert to a Lymphoma Phenotype

Abstract

The MYC transcription factor regulates a vast number of genes and is implicated in many human malignancies. In some hematological malignancies, MYC is frequently subject to missense mutations that enhance its transformation activity. Here, we use a novel murine cell system to (i) characterize the transcriptional effects of progressively increasing MYC levels as normal primary B-cells transform to lymphoma cells and (ii) determine how this gene regulation program is modified by lymphoma-associated MYC mutations (T58A and T58I) that enhance its transformation activity. Unlike many previous studies, the cell system exploits primary B-cells that are transduced to allow regulated MYC expression under circumstances where apoptosis and senescence pathways are abrogated by the over-expression of the Bcl-xL and BMI1 proteins. In such cells, transition from a normal to a lymphoma phenotype is directly dependent on the MYC expression level, without a requirement for secondary events that are normally required during MYC-driven oncogenic transformation. A generalized linear model approach allowed an integrated analysis of RNA sequencing data to identify regulated genes in relation to both progressively increasing MYC level and wild type or mutant status. Using this design, a total of 7569 regulated genes were identified, of which the majority (n = 7263) were regulated in response to progressively increased levels of wild type MYC, while a smaller number of genes (n = 917) were differentially regulated, compared to wild type MYC, in T58A MYC- and/or T58I MYC-expressing cells. Unlike most genes that are similarly regulated by both wild type and mutant MYC genes, the set of 917 genes did not significantly overlap with known lipopolysaccharide regulated genes, which represent genes regulated by MYC in normal B cells. The genes that were differently regulated in cells expressing mutant MYC proteins were significantly enriched in DNA replication and G2 phase to mitosis transition genes. Thus, mutants affecting MYC proteins may augment quantitative oncogenic effects on the expression of normal MYC-target genes with qualitative oncogenic effects, by which sets of cell cycle genes are abnormally targeted by MYC as B cells transition into lymphoma cells. The T58A and T58I mutations augment MYC-driven transformation by distinct mechanisms.

Keywords: DNA replication; MYC; MYC mutants (T58A and T58I); cell cycle; lymphoma; oncogene deregulation; ribosome biogenesis.

Figure 1

Phenotypic characterization of immortal B-cells expressing doxycycline regulatable MYC. (A) Experimental setup, time points for doxycycline wash out and re-addition as well as time points for data collection and RNA sampling (B) MYC protein levels at different levels of doxycycline in culture medium (0, 25, 50, 100, 300, 600 and 1000 ng/mL) assessed by Western Blot for WT MYC (red), T58A MYC (green) and T58I MYC (blue) in relation to actin levels. Mean ± SD (n = 3) is shown. (C) Flow cytometry data for WT MYC (red), T58A (green) and T58I (blue) sampled at 48, 96 and 144 h following addition of doxycycline at 0, 25, 50, 100, 300, 600 and 1000 ng/mL. From top to bottom, the panels represent: cumulative live cell number, relative cell size estimated by the forward scatter channel (FSC) and G₀/G₁, S and G₂/M cell cycle stage assessment for cells stained with propidium iodide (percentages). The data are from a single experiment. Similar data for the 48 h time point (used for RNA sequencing experiments) were obtained in a separate experiment.

Figure 2

Principal component analysis of global gene expression data. The first three principal components describe 75.6% of the variation in the global gene expression data across the three MYC genotypes and seven expression levels. Levels of MYC WT (circles), MYC T58A (triangles) and MYC T58I (squares) were controlled by titrating in doxycycline at the indicated concentrations (ng/mL, colors). (A) Scatter plot of principal component (PC) 1 versus PC2. (B) Scatter plot of PC2 versus PC3. The proportion (%) of the variation accounted for by each PC is shown.

Figure 3

Expression levels for gene clusters. Hierarchical clustering based on expression levels for the union of genes that were either significantly changed upon increased MYC levels in WT or that were differentially regulated in response to increasing MYC levels between MYC WT and MYC T58A or between MYC WT and MYC T58I (7569 genes) yielded 12 clusters (n = 3 independent experiments). The heatmap shows reads per kilobase per million of reads (RPKM) for all clustered genes (rows) for different levels (plotted above the heat map) of the three MYC genotypes (MYC WT (red dots), MYC T58A (green) and MYC T58I (blue), columns). Doxycycline concentrations (ng/mL) used to treat the cells are indicated by the bar-chart above the main heatmap. To the right of the main heatmap are three tracks describing, from left to right, whether the differentially regulated genes encode for protein-coding (yellow) or non-coding (light green) RNA, whether (blue) or not (pink) a gene has a canonical E-box motif within 1kb of its transcription start site (TSS), and the last track describes the distance between the TSS and the closest canonical E-box motif.

Figure 4

Gene ontology functional classification of clustered genes. Over-representation tests of clustered genes in gene sets with different gene ontology (GO) terms, representing biological processes, for clusters c1 and c3–c12. The most enriched gene sets following GO term simplification are shown. Median RPKM values for genes in denoted clusters are plotted (right panels) for MYC WT (red), MYC T58A (green) and MYC T58I (blue) for increasing doxycycline levels (0, 25, 50, 100, 300, 600 and 1000 ng/mL), and the colored lines corresponding to the respective MYC genotype represent a 2nd degree polynomial regression model where the shaded areas indicate the 95% confidence interval for the model (larger high resolution plots are available in Supplementary Figure S4).

Figure 5

Comparisons with published datasets relevant for MYC-driven lymphogenesis. (A) Overlap between differentially regulated genes in WT MYC and significantly changed genes between any of the conditions measured by Sabò et al. in the Eµ-Myc mouse model [18]. p-value from Fisher’s exact test. (B) Representation of significantly changed genes in the Sabò et al. dataset in gene clusters c1–c12. Genes from clusters c3, c7, c9 and c10 were over-represented in the regulated genes from the Sabò et al. study (Fisher’s exact test, p = 0.02, 6.0 × 10⁻²⁰, 2.0 × 10⁻²⁰ and 2.0 × 10⁻⁹, respectively). (C) Overlap between 88 of the 100 most significant direct MYC target genes derived from Muhar et al. [21] that could be converted to the mouse ortholog with a unique gene identifier and genes with significantly altered transcript levels in WT MYC. p-value from Fisher’s exact test. (D) A total of 6129 genes that were differentially regulated in WT MYC and that had unique gene identifiers following ortholog conversion between human and mouse ranked by adjusted p-value. Each grey bar represents one gene with the 71 genes that were also present in the data from Muhar et al., marked in black. The set of 71 genes was significantly enriched in the set of most significantly regulated genes in WT MYC (Fisher’s exact test, p-value = 0.0012).

Figure 6

Differences in differentially regulated genes between MYC WT, MYC T58A and MYC T58I. (A) Venn diagram representation of genes with significant transcript level changes (FDR q-value ≤ 0.01 and max fold change ≥ 2). The diagram shows the set of genes with significant regulation in MYC WT (WT) and its overlap with the sets of genes that are regulated significantly differently from WT MYC in MYC T58A (T58A vs. WT) and/or MYC T58I (T58I vs. WT). (B) Heatmap representation of Fisher’s exact test significance levels for intersects between the Venn diagram subsets in panel A (A–G) and genes in clusters 1–12. Intersects with a significantly larger overlap than expected are represented in red (over-rep) and intersects with an overlap significantly smaller than expected are represented in blue (under-rep). (C) Comparison differentially regulated gene sets (A and B–G) with lipopolysaccharide (LPS) regulated genes, representing normal MYC target genes. The names of LPS-regulated genes were obtained from previously published work [17]; supplementary data file “embr201947987-sup-0004-datasetev2.xlsx”, all groups. Venn diagrams show a significant overlap between Set A genes and LPS-regulated genes and a smaller overlap for Set B–G genes (all sets except for Set A) that is not statistically significant. p-values to test over-representation of genes in the overlaps are shown (Fisher’s exact test, one-tailed).

Figure 7

Enriched pathways for genes with sensitivity differences to the levels of MYC T58A and MYC T58I. Representative examples of gene sets that are significantly enriched in a gene set enrichment analysis (GSEA) of genes ranked according to the extent of their sensitivity to T58A MYC in relation to T58I MYC. Left panels: GSEA plots showing a running sum statistic that increases when genes in the ranked list are present in the functional gene set tested. The position of these genes in the rank is indicated by vertical lines. The red dashed line shows the enrichment score value (ES, maximum value of the running sum score). The analysis was performed on 11,737 genes for which human orthologues could reliably be identified and the sensitivity-based rank is shown in the lowest plot. The pathway name is indicated above each plot followed by the normalized enrichment score (NES), GSEA p-value, false discovery rate adjusted p-value (FDR) and the number of genes in the GSEA leading edge (interval between the start of the rank and the gene where ES is attained) in relation to the number of genes in the pathway. Middle panel: Expression-level plots showing median reads per kilobase per millions of reads (RPKM) values for leading edge set of genes identified for each GSEA gene set. Median RPKM values for WT MYC (red), T58A MYC (green) and T58I MYC (blue) at different doxycycline levels (0, 25, 50, 100, 300, 600 and 1000 ng/mL) are plotted. Right panels: Violin plots showing the MYC-level sensitivity of the leading-edge genes in the indicated pathway for WT MYC (red), T58A MYC (green) and T58I MYC (blue). Sensitivity was determined by the change in absolute fold change between the fold change between 0 and 25 ng/mL and the fold change between 300 and 600 ng/mL of doxycycline (FC_25/0/FC_600/300).

Figure 8

Genes that are differentially regulated by T58I MYC. (A) Plots showing the expression levels (RPKM) as a function of doxycycline concentration (0, 25, 50, 100, 300, 600 and 1000 ng/mL) for WT MYC (red), T58A MYC (green) and T58I MYC (blue). (A) Representative genes that show lower regulation by T58I MYC than WT MYC and/or T58A MYC. (B) Representative genes that show higher regulation by T58I MYC than WT MYC and/or T58A MYC. (C) Representative genes that show a different regulation direction by T58I MYC than WT MYC and/or T58A MYC. (D) Protein conformation predictions for the conserved MYC-Box I region, containing residue 58 (dashed line), in WT MYC (red), T58A MYC (green) and T58I MYC (blue). Dynamine predicts rigidity of the protein backbone (scale 0–1, where 1 is highly rigid). Anchor predicts protein interaction propensity (scale 0–1, where 1 is the highest propensity level). Helix (Agadir) is a prediction of alpha-helicity using the Agadir algorithm (the results are essentially identical for the three proteins, explaining the existence of only one visible plotted line).