Immunomodulatory drugs (IMiD) are a backbone therapy for multiple myeloma (MM). Despite their efficacy, most patients develop resistance, and the mechanisms are not fully defined. Here, we show that IMiD responses are directed by IMiD-dependent degradation of IKZF1 and IKZF3 that bind to enhancers necessary to sustain the expression of MYC and other myeloma oncogenes. IMiD treatment universally depleted chromatin-bound IKZF1, but eviction of P300 and BRD4 coactivators only occurred in IMiD-sensitive cells. IKZF1-bound enhancers overlapped other transcription factor binding motifs, including ETV4. Chromatin immunoprecipitation sequencing showed that ETV4 bound to the same enhancers as IKZF1, and ETV4 CRISPR/Cas9-mediated ablation resulted in sensitization of IMiD-resistant MM. ETV4 expression is associated with IMiD resistance in cell lines, poor prognosis in patients, and is upregulated at relapse. These data indicate that ETV4 alleviates IKZF1 and IKZF3 dependency in MM by maintaining oncogenic enhancer activity and identify transcriptional plasticity as a previously unrecognized mechanism of IMiD resistance.

Significance:

We show that IKZF1-bound enhancers are critical for IMiD efficacy and that the factor ETV4 can bind the same enhancers and substitute for IKZF1 and mediate IMiD resistance by maintaining MYC and other oncogenes. These data implicate transcription factor redundancy as a previously unrecognized mode of IMiD resistance in MM.

See related article by Welsh, Barwick, et al., p. 34.

See related commentary by Yun and Cleveland, p. 5.

This article is featured in Selected Articles from This Issue, p. 4

Over the past two decades, the outcomes of patients with plasma cell malignancy multiple myeloma (MM) have dramatically improved because of the development of novel agents, such as immunomodulatory drugs (IMiD). IMiDs include thalidomide and its analogues lenalidomide, pomalidomide, and others (1, 2). IMiDs are reported to modulate the immune system, reduce angiogenesis, and cause cytostatic inhibition of MM cells (3); the mechanism underlying these activities was not discovered until after their clinical efficacy was established (4, 5). The IMiD thalidomide was first shown to bind the E3 ubiquitin ligase Cereblon (CRBN), which is necessary for its infamous teratogenic effects in utero (6). Subsequently, CRBN was shown to be necessary for the antimyeloma effects of lenalidomide (7), and this was due, at least in part, to IMiD-directed CRBN binding and subsequent ubiquitination and proteasomal degradation of the transcription factors IKAROS (IKZF1) and AIOLOS (IKZF3; refs. 8, 9). Before IKZF1 and IKZF3 were known to be targeted by IMiDs, Ikzf1 was identified as essential for lymphocyte commitment in the mouse (10). Later homology screens identified Ikzf3, which was subsequently shown to be necessary for plasma cell differentiation (11, 12). However, how IKZF1 and IKZF3, and therefore IMiDs, regulate the gene-expression program in MM cells remains poorly understood.

Gene regulation has been extensively studied in B cells, where immunoglobulin heavy (IGH) enhancers were the first eukaryotic enhancers discovered (13–15) and the gene regulatory program of plasma cells has been well characterized (16). Recent work has identified clustered enhancer elements that form super-enhancers (SE), which are densely bound by transcriptional coactivators such as the acetyltransferase P300, the RNA polymerase II elongation factor BRD4, and mediator of RNA polymerase II transcription subunit 1 (MED1; ref. 17). SEs control cell identity in normal and cancerous tissues, including MM cell lines (18) and the gene regulatory landscape of primary MM is characterized by SEs that drive oncogene expression (19–21). IKZF1 and IKZF3 have been shown to function at SEs in B-cell leukemia (22), but there are little data on the genomic targets of IKZF1 and IKZF3 in MM or on how IMiDs may regulate SEs.

Structural rearrangements in MM juxtapose SEs with oncogenes. Primary events include translocation of the IGH enhancers (23) to one of the Cyclin D genes (CCND1, CCND2, CCND3; ref. 24), the histone methyltransferase NSD2 (also known as MMSET and WHSC1), which often involves reciprocal translocation of the IGH μ enhancer to drive ectopic expression of FGFR3 (25, 26), or a Maf proto-oncogene (MAF, MAFB; ref. 27). Although almost always a secondary event, MYC is the most common dysregulated oncogene in MM and occurs by focal amplification of MYC or its enhancers and complex structural rearrangements, often involving IGH enhancers, immunoglobulin light chain (IGL, IGK) enhancers, SEs at TXNDC5, FAM46C, or other loci (28–33). MYC transcriptional amplification is associated with disease progression (32, 34) and is observed in most MM cell lines (35). IMiD sensitivity corresponds with the downregulation of MYC and IRF4, suppressing proliferation and cell survival, respectively (36, 37), but whether IKZF1 and IKZF3 directly regulate MYC is still not known.

In the present study, we identified IKZF1 binding to the majority of enhancers and almost all SEs throughout the MM epigenome. Although IMiDs universally depleted chromatin-bound IKZF1, the transcriptional coactivators BRD4 and P300 were selectively depleted in IMiD-sensitive cells, but not in IMiD-resistant cells. IMiD-resistant MM cells maintained BRD4 and P300 occupancy and oncogenic enhancer function through the expression of the E26 transformation-specific (ETS) transcription factor ETV4, which cobinds enhancers with IKZF1 and induces IMiD resistance. These results demonstrate that transcriptional plasticity is a previously unappreciated mechanism of IMiD resistance in MM.

IMiDs Downregulate MYC in IMiD-Sensitive MM Cells

IMiD activity was measured in 12 MM cell lines treated with lenalidomide for up to six days (Fig. 1A). This analysis identified IMiD-sensitive (e.g., MM1S, OPM2, and H929) and resistant cells (e.g., RPMI8226, L363, KMS18, INA6, and JJN3). Because IMiD sensitivity is reported to correspond with IRF4 and MYC downregulation (36, 37), mRNA expression was determined by RT-qPCR and indicated that IRF4 and MYC were repressed in lenalidomide-sensitive cells, but not in lenalidomide-resistant cells, where compensatory transcriptional upregulation often occurred (Supplementary Fig. S1A). MYC protein depletion was also confirmed in sensitive MM cells but not in resistant cells, despite IKZF1 and IKZF3 depletion (Fig. 1B; Supplementary Fig. S1B). As expected, lenalidomide-mediated depletion of IKZF3 and MYC was dependent on CRBN, as demonstrated in CRBN-deficient MM cells (Supplementary Fig. S1C), which is consistent with previous results (7).

Figure 1.

Lenalidomide (Len) downregulates MYC in Len-sensitive MM cells. A, MTT assay of MM cells treated with Len (10 μmol/L) at the indicated time points. The line shows the threshold separating sensitive (Len sen) from resistant (Len res) cell lines. B, Immunoblot analysis of IKZF3, MYC, CRBN, and GAPDH in 7 MM cell lines treated with lenalidomide (10 μmol/L) at the indicated time points grouped by lenalidomide response. C, Scatter plot of RNA-seq gene-expression changes in MM1S (top) and OPM2 (bottom) treated with lenalidomide for 12 and 24 hours, respectively. Differentially expressed genes (DEG; FDR ≤0.05) are shown in red (upregulated) and blue (downregulated). Ctrl, control; FPKM, Fragments Per Kilobase Per Million reads. D, Heat map of lenalidomide DEGs in MM1S (top) and OPM2. Samples are represented in columns and genes by rows. E, Gene set enrichment analysis (GSEA) enrichment scores for the gene sets: HALLMARK_INTERFERON_ALPHA_RESPONSE (top) and MYC_TARGET_V1 (bottom) in both MM1S (blue) and OPM2 (green). Lenalidomide-induced gene-expression changes are sorted from most upregulated (left) to most downregulated (right) and the overlap with gene set genes are shown in blue and green for MM1S and OPM2, respectively. F, Expression of selected DEGs measured by RNA-seq. Data presented are the mean ± SEM (A) or ± SD (F) with N ≥ 3 (A) N = 2 (F).

Figure 1.

Lenalidomide (Len) downregulates MYC in Len-sensitive MM cells. A, MTT assay of MM cells treated with Len (10 μmol/L) at the indicated time points. The line shows the threshold separating sensitive (Len sen) from resistant (Len res) cell lines. B, Immunoblot analysis of IKZF3, MYC, CRBN, and GAPDH in 7 MM cell lines treated with lenalidomide (10 μmol/L) at the indicated time points grouped by lenalidomide response. C, Scatter plot of RNA-seq gene-expression changes in MM1S (top) and OPM2 (bottom) treated with lenalidomide for 12 and 24 hours, respectively. Differentially expressed genes (DEG; FDR ≤0.05) are shown in red (upregulated) and blue (downregulated). Ctrl, control; FPKM, Fragments Per Kilobase Per Million reads. D, Heat map of lenalidomide DEGs in MM1S (top) and OPM2. Samples are represented in columns and genes by rows. E, Gene set enrichment analysis (GSEA) enrichment scores for the gene sets: HALLMARK_INTERFERON_ALPHA_RESPONSE (top) and MYC_TARGET_V1 (bottom) in both MM1S (blue) and OPM2 (green). Lenalidomide-induced gene-expression changes are sorted from most upregulated (left) to most downregulated (right) and the overlap with gene set genes are shown in blue and green for MM1S and OPM2, respectively. F, Expression of selected DEGs measured by RNA-seq. Data presented are the mean ± SEM (A) or ± SD (F) with N ≥ 3 (A) N = 2 (F).

Close modal

To more broadly characterize the transcriptional effects of IMiDs in MM, RNA sequencing (RNA-seq) was performed on IMiD-sensitive MM1S and OPM2 cells treated with lenalidomide. This analysis identified 320 and 1,239 differentially expressed genes (DEG) in MM1S and OPM2, respectively (Fig. 1C; Supplementary Data 1). DEGs were consistent between replicates in each cell type, as indicated by hierarchical clustering, which separated the lenalidomide-treated from mock-treated controls (Fig. 1D). Annotation of gene-expression changes by gene set enrichment analysis (GSEA) indicated that lenalidomide increased interferon-α gene expression in both MM1S and OPM2 (Fig. 1E, top), consistent with previous reports (36, 38). Upregulated interferon-α genes included IFIT3 and MX2, as well as other type I interferon responsive genes not included in the GSEA Hallmark Interferon-α gene set, such as CD38 and CD274 (which encodes PD-L1; Fig. 1F, top; refs. 36, 39). Downregulated genes included MAF in MM1S, which has an IGH;MAF t(14;16) translocation, and FGFR3 in OPM2, which has an IGH;FGFR3 t(4;14) translocation (Fig. 1F, middle). Together, these data suggest that IMiDs impair the function of the IGH enhancer(s), which are translocated to these oncogenes. This was further supported by the downregulation of MYC in MM1S, which also has an IGH;MYC t(8;14) translocation (Fig. 1F, bottom). Furthermore, lenalidomide also resulted in MYC downregulation in OPM2 (Fig. 1F, bottom), and MYC targets were the most depleted gene set of all Hallmark gene sets in both cell types (Fig. 1E, bottom; Supplementary Data S1).

The above results were corroborated by silencing IKZF3 using shRNAs. Two separate OPM2 shIKZF3 clones were established, and doxycycline was used to induce IKZF3 knockdown. RNA-seq analysis of doxycycline-treated cells identified 298 dysregulated genes, including IKZF3 (Supplementary Fig. S2A; Supplementary Data S2). The results were similar between both OPM2 clones, as demonstrated by the hierarchical clustering of DEGs (Supplementary Fig. S2B). Similar to the gene-expression analysis of lenalidomide-treated cells, GSEA indicated upregulation of type I interferon genes and downregulation of MYC signaling (Supplementary Fig. S2C; Supplementary Data S2). Furthermore, GSEA of lenalidomide-regulated genes in OPM2 (Fig. 1C) showed that lenalidomide-upregulated genes were enriched in the shRNA controls and that lenalidomide-downregulated genes were enriched in the shIKZF3 OPM2 cells (Supplementary Fig. S2D). Taken together, these results indicate that IMiDs downregulate MYC by depleting IKZF1 and/or IKZF3 and suggest that a “bottleneck” in IMiD resistance cells is MYC downregulation rather than IKZF1 and IKZF3 degradation.

IKZF1 Binds MYC and MM Super-Enhancers

IKZF1 chromatin immunoprecipitation sequencing (ChIP-seq) was performed on MM1S cells to identify direct targets of IMiDs. This analysis identified 45,293 IKZF1 binding sites relative to those in the input control (Fig. 2A; Supplementary Data S3). Inspection of genomic regions near lenalidomide-dysregulated genes identified IKZF1 binding at the MYC locus, including both regions immediately proximal to MYC (Supplementary Fig. S3A), as well as a cluster of regions ∼600 kb distal to MYC on the telomeric side of 8q24, which coincides with a cluster of translocations and focal copy-number gains found in primary MM samples from the CoMMpass study (Fig. 2B). Although other genes were more proximal to this translocation hotspot, MYC was the only gene downregulated by lenalidomide treatment in OPM2 and MM1S, consistent with recent data that suggest that enhancers in this region regulate MYC (BioRxiv 2023.05.19.541506). Notably, MM1S contains a t(8;14) MYC;IGH translocation telomeric of the region shown here. IKZF1 was also directly bound to several type I interferon-regulated genes, including IFIT3, CD274, and CD38 (Supplementary Fig. S3A). These data indicate that IKZF1-bound loci have both repressive and activating functions, as previously reported in leukemia (22).

Figure 2.

IKZF1 binds MM enhancers and super-enhancers. A, Heat map of IKZF1 of ChIP-seq peaks (left) and the input control (right) in MM1S cells. The scale is expressed in reads per million (RPM). B, IKZF1 binding at the MYC (red) locus compared with the input control. Translocation breakpoints (Tr, black) and the frequency of copy-number variant gains (CNV; gray) from CoMMpass samples are shown (bottom). C, P300 (blue) super-enhancer analysis in MM1S (top) and RPMI8226 (bottom). Binding of IKZF1 (orange), H3K27ac (gray), and BRD4 (green) to each P300-defined region is shown using the normalized RPM. D, Heat maps of P300, H3K27ac, BRD4, and IKZF1 in P300-defined enhancers from MM1S (top) and RPMI8226 (bottom). All heat maps were sorted by the size of the P300-defined super-enhancer (SE; annotated left). E, Overlap of H3K27ac, BRD4, and IKZF1 with P300 enhancers (blue) and super-enhancers (black). F, Genome plots of CCND2 (left) and BCL2L1 (right) loci showing P300 (blue), H3K27ac (gray), BRD4 (green), and IKZF1 (orange) binding in MM1S (top) and RPMI8226 (bottom). Scale is RPM.

Figure 2.

IKZF1 binds MM enhancers and super-enhancers. A, Heat map of IKZF1 of ChIP-seq peaks (left) and the input control (right) in MM1S cells. The scale is expressed in reads per million (RPM). B, IKZF1 binding at the MYC (red) locus compared with the input control. Translocation breakpoints (Tr, black) and the frequency of copy-number variant gains (CNV; gray) from CoMMpass samples are shown (bottom). C, P300 (blue) super-enhancer analysis in MM1S (top) and RPMI8226 (bottom). Binding of IKZF1 (orange), H3K27ac (gray), and BRD4 (green) to each P300-defined region is shown using the normalized RPM. D, Heat maps of P300, H3K27ac, BRD4, and IKZF1 in P300-defined enhancers from MM1S (top) and RPMI8226 (bottom). All heat maps were sorted by the size of the P300-defined super-enhancer (SE; annotated left). E, Overlap of H3K27ac, BRD4, and IKZF1 with P300 enhancers (blue) and super-enhancers (black). F, Genome plots of CCND2 (left) and BCL2L1 (right) loci showing P300 (blue), H3K27ac (gray), BRD4 (green), and IKZF1 (orange) binding in MM1S (top) and RPMI8226 (bottom). Scale is RPM.

Close modal

Given that downregulation of MYC and other essential genes correlated with lenalidomide responses, we investigated the IKZF1-bound loci associated with gene activation. To do this, ChIP-seq was performed for the transcriptional coactivators P300 and BRD4, as well as the enhancer histone modification H3K27ac, in IMiD-sensitive MM1S cells. A comparison of the aforementioned enhancer marks indicated that 51% of IKZF1-bound regions overlapped with P300 binding sites (Supplementary Fig. S3B, left), implying that approximately half of the IKZF1 sites were enhancers. Because IKZF1 is reported to regulate SEs in leukemia (22), we determined whether IKZF1 was bound to such regions in MM using P300 to stitch together enhancers and define SEs. This analysis identified 687 P300-defined SEs in MM1S, including those previously reported in MM, such as DUSP22/IRF4, IGH, IGL, CCND2, MCL1, PRDM1, BCL2L1 (which encodes BCL-XL); structural variant MM hotspots (29, 30) including IGH, MYC, IGL, and TXNDC5; and MM dependencies such as IRF4, CCND2, MEF2C, IRF2, POU2AF1 (Fig. 2C, top; Supplementary Data S3; ref. 40). Overlaying of IKZF1 binding levels in these regions suggested that IKZF1 was cobound at the P300 sites, including the aforementioned SEs (Fig. 2C, see orange). This was evident by comparing a pile-up of ranked P300 regions with H3K27ac, BRD4, and IKZF1 binding, which indicated that the majority of the regions were common (Fig. 2D, top). Indeed, the quantification of P300 and IKZF1 overlap showed that IKZF1 was found at ∼80% of the P300 sites and virtually at all P300-defined SEs (Fig. 2E, top). This included SEs at CCND2 and BCL2L1, demonstrating that P300, BRD4, and IKZF1 occupied the same genetic elements (Fig. 2F, top).

To determine whether IKZF1 binding to enhancers and SEs was specific to IMiD-sensitive MM cells, lenalidomide-resistant RPMI8226 cells were also characterized by ChIP-seq for IKZF1, P300, BRD4, and H3K27ac. The number of P300-defined enhancers (N = 14,136) and SEs (N = 793) in RPMI8226 were similar to those found in MM1S, and the P300 sites tended to be cobound with H3K27ac, BRD4, and IKZF1, similar to MM1S (Fig. 2C and D, bottom; Supplementary Fig. S3B, right). RPMI8226 IKZF1 bound elements overlapped ∼80% of the P300 enhancers and nearly all P300-defined SEs, including those at the CCND2 and BCL2L1 loci (Fig. 2E and F, bottom). Because IMiD deplete both IKZF1 and IKZF3 (8, 9), we also investigated IKZF3 binding using cleavage under targets and tagmentation (CUT&Tag; ref. 41) for IKZF1, IKZF3, and the histone modifications H3K27ac, H3K4me3, and H3K27me3. Clustering of IKZF1- or IKZF3-bound regions in both cell lines (MM1S and RPMI8226) corroborated the IKZF1 ChIP-seq data showing that approximately half of these sites overlapped H3K27ac regions (Supplementary Fig. S4A), including those at the IRF4 and IGL SEs (Supplementary Fig. S4B). These data indicate that IKZF1 and IKZF3 are bound to enhancers and SEs in MM regardless of lenalidomide sensitivity.

Lenalidomide Displaces Transcriptional Cofactors in IMiD-Sensitive MM Cells

Given that the majority of enhancers were bound by IKZF1 in both IMiD-sensitive and resistant MM cells, we next sought to determine whether lenalidomide differentially affected these enhancers. MM1S and RPMI8226 cells were treated with 10 μmol/L lenalidomide prior to performing IKZF1 ChIP-seq analysis. As expected, in IMiD-sensitive MM1S cells, lenalidomide reduced the amount of IKZF1-bound chromatin (Fig. 3A). To further determine how lenalidomide treatment affected coactivator function, P300 and BRD4 ChIP-seq were also performed in MM1S cells treated with lenalidomide. As shown in Fig. 3B, P300 was dysregulated compared with that in the control and overall depleted. Notably, BRD4 was more significantly reduced (Fig. 3C), suggesting that the loss of IKZF1 resulted in reduced enhancer activity. In IMiD-resistant RPMI8226 cells, lenalidomide also resulted in a significant reduction in chromatin-bound IKZF1 (Fig. 3D), indicating that lenalidomide resistance in RPMI8226 cells was not due to failure to deplete IKZF1, consistent with western blot analysis (Supplementary Fig. S1B). However, inspection of P300 and BRD4 binding in lenalidomide-treated RPMI8226 showed that P300 and BRD4 were very similar to control cells and not significantly reduced (Fig. 3E and F), suggesting that the enhancer function was intact despite IKZF1 depletion.

Figure 3.

Lenalidomide-induced loss of enhancer coactivators in IMiD-sensitive cells. AC, Scatter plots of IKZF1 (A), P300 (B), and BRD4 (C) ChIP-seq signal in MM1S control (ctrl) and lenalidomide (Len)-treated cells (10 μmol/L, 24 hours). DF, Scatter plots of IKZF1 (D), P300 (E), and BRD4 (F) in RPMI8226 ctrl and Len-treated cells. The signal at each region is shown to the right of the scatter plots for control and lenalidomide treatment. G, Scatter plot of lenalidomide-induced fold changes at regions cobound with IKZF1 and P300 (left), as well as a binned analysis of fold changes ranked by lenalidomide-induce IKZF1 fold change (right) in MM1S. Correlation is shown (Spearman ρ). H, Scatter plot and binned analysis of IKZF1 and BRD4 cobound regions for MM1S as in part G. I, Scatter plot and binned analysis of IKZF1 and P300 cobound regions for RPMI8226 as in part G. J, Scatter plot and binned analysis of IKZF1 and BRD4 cobound regions for RPMI8226 as in G. K, Genome plot of DUSP22 and IRF4 showing IKZF1 (orange), P300 (blue), and BRD4 (green) ChIP-seq in MM1S (left) and RPMI8226 (right). The change in lenalidomide-treated samples is shown in darker tones, and the scale is reads per million (RPM).

Figure 3.

Lenalidomide-induced loss of enhancer coactivators in IMiD-sensitive cells. AC, Scatter plots of IKZF1 (A), P300 (B), and BRD4 (C) ChIP-seq signal in MM1S control (ctrl) and lenalidomide (Len)-treated cells (10 μmol/L, 24 hours). DF, Scatter plots of IKZF1 (D), P300 (E), and BRD4 (F) in RPMI8226 ctrl and Len-treated cells. The signal at each region is shown to the right of the scatter plots for control and lenalidomide treatment. G, Scatter plot of lenalidomide-induced fold changes at regions cobound with IKZF1 and P300 (left), as well as a binned analysis of fold changes ranked by lenalidomide-induce IKZF1 fold change (right) in MM1S. Correlation is shown (Spearman ρ). H, Scatter plot and binned analysis of IKZF1 and BRD4 cobound regions for MM1S as in part G. I, Scatter plot and binned analysis of IKZF1 and P300 cobound regions for RPMI8226 as in part G. J, Scatter plot and binned analysis of IKZF1 and BRD4 cobound regions for RPMI8226 as in G. K, Genome plot of DUSP22 and IRF4 showing IKZF1 (orange), P300 (blue), and BRD4 (green) ChIP-seq in MM1S (left) and RPMI8226 (right). The change in lenalidomide-treated samples is shown in darker tones, and the scale is reads per million (RPM).

Close modal

We next compared the IMiD-mediated IKZF1 and P300 changes at sites cobound by both factors in MM1S cells and observed a strong correlation (Spearman ρ = 0.481; Fig. 3G, left), indicating that the more IKZF1 was depleted, the more chromatin-bound P300 was lost. This could also be observed by binning the IMiD-induced changes in P300 binding according to the fold change of IKZF1 (Fig. 3G, right), which showed that the amount of IKZF1 loss correlated with the amount of P300 loss at the cobound sites. A similar analysis performed with the BRD4 and IKZF1 cobound sites revealed the same phenomenon (Fig. 3H). In contrast, there was no correlation with P300 or BRD4 changes in IMiD-resistant RPMI8226 cells (Fig. 3I and J). This was also observed at the DUSP22/IRF4 SE, where IMiD-induced loss of IKZF1, P300, and BRD4 in MM1S cells, as depicted by the negative values for the change in binding upon lenalidomide treatment (Fig. 3K, left). Similarly, in RPMI8226, lenalidomide resulted in reduced IKZF1 binding. However, P300 and BRD4 did not exhibit reduced binding compared with MM1S (Fig. 3K, right). This was also true for the IGH enhancer, which drives MAF and MYC expression in MM1S and showed lenalidomide-induced loss of P300 and BRD4 in MM1S but not in RPMI8226 (Supplementary Fig. S5). Taken together, these data demonstrate that lenalidomide treatment is associated with release of transcriptional coactivators from chromatin in lenalidomide-sensitive MM1S cells, but this fails to occur in lenalidomide-resistant RPMI8226 cells.

ETV4 Binds IKZF1 Enhancers in Lenalidomide-Resistant MM Cells

Because transcription factors can mediate the recruitment of coactivators, including P300 and BRD4 (42), we investigated whether other transcription factors could be compensating for the loss of IKZF1 upon lenalidomide treatment in RPMI8226 cells. First, we determined the binding motifs enriched in the IKZF1-bound regions and identified 135 motifs (Supplementary Data S4; FDR ≤10−10, odds ratio ≥1.5, RNA FPKM ≥1). These transcription factor binding motifs were clustered based on the overlap with specific IKZF1-bound regions, such that motifs cooccurring at the same subsets of IKZF1 sites were highly correlated with each other. This revealed families of transcription factor motifs present in distinct sets of IKZF1-bound elements, including multiple classes of zinc finger motifs, ETS, immunoglobulin enhancer-binding factors E12 and E47 (E2A), basic helix–loop–helix leucine zipper (bHLH-ZIP), interferon regulatory factor (IRF), and activator protein-1 (AP-1) motifs (Fig. 4A, left). The fraction of IKZF1 sites overlapped by each motif ranged from 1.0% to 70.9%, with the highest averages observed in the E2A (28.2%), ETS (27.1%), zinc finger (18.3%), bHLH-ZIP (15.3%), activator protein-1 (AP-1; 11.6%), and IRF (4.0%) families (Fig. 4A, right). Interestingly, the canonical IKZF1 motif clustered with ETS factors, away from other zinc finger factors (Fig. 4A). Notably, we observed that IKZF1 and ETS factors (e.g., ETV4) shared a common “AGGAA” motif, which was not present in the other motif families (Fig. 4B). We further prioritized transcription factors by correlating basal RNA expression with response to lenalidomide in the 12 MM cell lines characterized in Fig. 1A as well as with progression-free survival (PFS) for patients treated with IMiD-containing regimens on the Multiple Myeloma Research Foundation (MMRF) CoMMpass trial. This analysis identified several factors correlating with lenalidomide resistance in cell lines, including the Kruppel-like factors KLF4 and KLF5, which contain C2H2 zinc finger domains; AP-1 factors JUNB and FOS; as well as the ETS factor ETV4 (Fig. 4C; Supplementary Data S4). We chose to focus on ETV4, given its binding motif similarity to IKZF1 (Fig. 4B), its association with poor PFS in primary samples (Fig. 4D), and its high correlation with lenalidomide resistance in MM cell lines (Fig. 4E). To experimentally determine whether ETV4 and IKZF1 bound to the same elements, we performed ChIP-seq for ETV4 in RPMI8226 and identified 15,678 ETV4-enriched regions relative to the input control (Fig. 5A; Supplementary Data S5). By comparing IKZF1 and ETV4 genomic binding sites, we observed that the occupancy of these two factors was highly correlated (Fig. 5B) with 78.9% (N = 12,364) of ETV4 sites overlapping IKZF1-bound regions (Fig. 5C). Interestingly, the IKZF1 sites that overlapped ETV4 regions had a higher prevalence of the cofactors P300 and BRD4 as well as the enhancer modification H3K27ac compared with other IKZF1 sites (Fig. 5D, compared with Supplementary Fig. S3B). Furthermore, an inspection of SE loci at IGH and PIM1 indicated that ETV4 bound the same genomic elements as IKZF1, BRD4, and P300 (Fig. 5E). Taken together, these data suggest that ETV4, when expressed in MM, functions at IKZF1-bound elements as an extralineage transcription factor sustaining oncogenic enhancers.

Figure 4.

IKZF1 shares similar motifs with ETS factors that correlate with IMiD resistance. A, Correlation matrix of motifs at IKZF1-bound regions in RPMI8226 (left) with select transcription factor families annotated in color (key far left). The frequency of overlap with IKZF1-bound regions is shown (right). Only motifs for expressed (≥1 FPKM) transcription factors in RPMI8226 significantly (FDR ≤10−10, odds ratio ≥1.5) enriched at IKZF1-bound regions are shown. B, Logo plots of select motifs enriched at IKZF1-bound regions in RPMI8226. The frequency of overlap with IKZF1-bound regions is shown in parenthesis. C, Pearson correlation (R) of transcription factor expression as measured by RNA-seq with lenalidomide (Len; 10 μmol/L)-induced change in growth index on Day 6 as measured by MTT assay in 12 MM cell lines from Fig. 1A. D, Progression-free survival (PFS) and hazard ratio (HR) for the expression of each transcription factor shown in C and observed in IMiD-treated patients from CoMMpass. E, Correlation of ETV4 expression as determined by RNA-seq with Len-induced change in growth index in 12 MM cell lines as determined by MTT assay on Day 6 (from Fig. 1A).

Figure 4.

IKZF1 shares similar motifs with ETS factors that correlate with IMiD resistance. A, Correlation matrix of motifs at IKZF1-bound regions in RPMI8226 (left) with select transcription factor families annotated in color (key far left). The frequency of overlap with IKZF1-bound regions is shown (right). Only motifs for expressed (≥1 FPKM) transcription factors in RPMI8226 significantly (FDR ≤10−10, odds ratio ≥1.5) enriched at IKZF1-bound regions are shown. B, Logo plots of select motifs enriched at IKZF1-bound regions in RPMI8226. The frequency of overlap with IKZF1-bound regions is shown in parenthesis. C, Pearson correlation (R) of transcription factor expression as measured by RNA-seq with lenalidomide (Len; 10 μmol/L)-induced change in growth index on Day 6 as measured by MTT assay in 12 MM cell lines from Fig. 1A. D, Progression-free survival (PFS) and hazard ratio (HR) for the expression of each transcription factor shown in C and observed in IMiD-treated patients from CoMMpass. E, Correlation of ETV4 expression as determined by RNA-seq with Len-induced change in growth index in 12 MM cell lines as determined by MTT assay on Day 6 (from Fig. 1A).

Close modal
Figure 5.

ETV4 is cobound at IKZF1 enhancers. A, Heat map of ETV4 ChIP-seq in RPMI8226 cells. B, Scatter plot of genomic binding occupancy for IKZF1 and ETV4 in RPMI8226 cells with color denoting if the region is bound by IKZF1 only, ETV4 only, or both (key right). C, Venn diagram of IKZF1 (orange) and ETV4 (red) bound regions in RPMI8226 cells. Note: the 11,847 ETV4-bound regions overlap 12,419 IKZF1-bound regions. D, Overlap of ETV4 sites with IKZF1, H3K27ac, P300, and BRD4 enriched regions in RPMI8226. E, Genomic plot of P300, H3K27ac, BRD4, IKZF1, and ETV4 at the IGH (left) and PIM1 (right) super-enhancers.

Figure 5.

ETV4 is cobound at IKZF1 enhancers. A, Heat map of ETV4 ChIP-seq in RPMI8226 cells. B, Scatter plot of genomic binding occupancy for IKZF1 and ETV4 in RPMI8226 cells with color denoting if the region is bound by IKZF1 only, ETV4 only, or both (key right). C, Venn diagram of IKZF1 (orange) and ETV4 (red) bound regions in RPMI8226 cells. Note: the 11,847 ETV4-bound regions overlap 12,419 IKZF1-bound regions. D, Overlap of ETV4 sites with IKZF1, H3K27ac, P300, and BRD4 enriched regions in RPMI8226. E, Genomic plot of P300, H3K27ac, BRD4, IKZF1, and ETV4 at the IGH (left) and PIM1 (right) super-enhancers.

Close modal

ETV4 Promotes IMiD Resistance

To further investigate the potential role of ETV4 in lenalidomide resistance, ETV4 ChIP-seq was conducted in lenalidomide-treated RPMI8226 cells. Similar to P300 and BRD4 in RPMI8226 cells, ETV4 genomic binding was largely unaffected by lenalidomide treatment (Supplementary Fig. S6A). Moreover, when we plotted the change in IKZF1 to the change in ETV4 at sites bound by both factors, we observed no correlation between IKZF1 and the loss of ETV4 (Supplementary Fig. S6B, left). This was further confirmed by plotting the change in ETV4 binned with the change in IKZF1 (Supplementary Fig. S6B, right). Notably, ETV4, P300, and BRD4 remained bound to the IGL enhancer following lenalidomide treatment despite IKZF1 depletion in lenalidomide-resistant RPMI8226 cells (Supplementary Fig. S6C, left), where the IGL enhancer is translocated with MAF and MYC oncogenes. In contrast, lenalidomide-sensitive MM1S cells lost IKZF1, P300, and BRD4 binding at the IGL enhancer (Supplementary Fig. S6C, right).

To determine whether ETV4 contributes to lenalidomide resistance in RPMI8226 cells, ETV4 was genetically disrupted using CRISPR/Cas9 genome editing. ETV4 protein levels were significantly reduced in RPMI8226 cells transduced with three different sgRNAs (ETV4 1–3) as compared with cells transduced with a nontargeting (NT) control (Fig. 6A). Reduced ETV4 expression was also observed at the mRNA level using RT-qPCR (Fig. 6B). Next, we treated RPMI8226 ETV4 CRISPR-edited cells with lenalidomide and observed that whereas IKZF3 was degraded in all lenalidomide-treated cells, MYC was downregulated only in ETV4-deficient cells and not in cells transduced with a nontargeted sgRNA (Fig. 6C). Consistent with IMiD-induced MYC downregulation, ETV4 knockout RPMI8226 cells were sensitized to lenalidomide whereas Cas9 control cells were not (Fig. 6D). We confirmed these data in two other lenalidomide-resistant cell lines that express ETV4: L363 and ARD. L363 has an IGL to MAP3K14 (encodes NIK) translocation, and following lenalidomide treatment, a pronounced downregulation of NIK was observed only in the ETV4KO (Fig. 6E), suggesting that the IGL enhancer is inhibited by lenalidomide only upon ETV4 depletion, as was the case in RPMI8226. As a result, the ETV4KO L363 became sensitive to lenalidomide (Fig. 6F). Similarly, the ETV4KO of ARD cells resulted in MYC downregulation upon lenalidomide treatment and were more sensitive to lenalidomide when compared with Cas9 control cells (Fig. 6G and H).

Figure 6.

ETV4 mediates lenalidomide resistance in MM cells. A, Western blot of ETV4 and GAPDH expression in MM1S, RPMI8226 with sgRNAs targeting ETV4 or a nontargeting (NT) control, and AMO1 MM cells. B,ETV4 RNA expression by RT-qPCR in RPMI8226 control cells (Cas9) and ETV4KO (sgETV4-1, sgETV-2, and sgETV4-3). C, Western blot for MYC, IKZF3, and GAPDH in RPMI8226 nontargeting control or sgRNA ETV4-2 and ETV4-3 cells treated with lenalidomide (Len; 10 μmol/L). D, F, H, Change in growth index measured by MTT assay at the indicated time points in RPMI8226 (D), L363 (F) and ARD (H) cell lines transduced with Cas9 control or sgRNA ETV4-2 and ETV4-3 cells after treatment with lenalidomide (Len; 10 μmol/L). E, G, Western blot analysis for NIK/MAP3K14, NFkB subunit (p100), ETV4, and MYC in L363 (E) and ARD cell lines (G) transduced with Cas9 control and ETV4KO cells treated with 10 μmol/L Len. Note the downregulation of NIK/MAP3K14 and lack of p100 degradation in ETV4KO cells in response to lenalidomide. F, Viability of L363 Cas9 and ETV4KO after Len treatment. The data presented are mean ± SEM.

Figure 6.

ETV4 mediates lenalidomide resistance in MM cells. A, Western blot of ETV4 and GAPDH expression in MM1S, RPMI8226 with sgRNAs targeting ETV4 or a nontargeting (NT) control, and AMO1 MM cells. B,ETV4 RNA expression by RT-qPCR in RPMI8226 control cells (Cas9) and ETV4KO (sgETV4-1, sgETV-2, and sgETV4-3). C, Western blot for MYC, IKZF3, and GAPDH in RPMI8226 nontargeting control or sgRNA ETV4-2 and ETV4-3 cells treated with lenalidomide (Len; 10 μmol/L). D, F, H, Change in growth index measured by MTT assay at the indicated time points in RPMI8226 (D), L363 (F) and ARD (H) cell lines transduced with Cas9 control or sgRNA ETV4-2 and ETV4-3 cells after treatment with lenalidomide (Len; 10 μmol/L). E, G, Western blot analysis for NIK/MAP3K14, NFkB subunit (p100), ETV4, and MYC in L363 (E) and ARD cell lines (G) transduced with Cas9 control and ETV4KO cells treated with 10 μmol/L Len. Note the downregulation of NIK/MAP3K14 and lack of p100 degradation in ETV4KO cells in response to lenalidomide. F, Viability of L363 Cas9 and ETV4KO after Len treatment. The data presented are mean ± SEM.

Close modal

To gain insight into how ETV4KO results in lenalidomide sensitization, we analyzed RPMI8226 and L363 Cas9 control and ETV4KO cells using the assay for transposase accessible chromatin-sequencing (ATAC-seq; ref. 43) in the presence or absence of lenalidomide treatment. Both ETV4KO and lenalidomide treatment resulted in consistent chromatin changes in both Cas9 and ETV4KO cells (Supplementary Fig. S7A and S7B). Chromatin accessibility gains with ETV4KO were consistently enriched for AP-1 motifs in both RPMI8226 and L363 (Supplementary Fig. S7C and S7D). Somewhat surprisingly, lenalidomide treatment primarily resulted in chromatin accessibility gains, which were consistent in both Cas9 control and ETV4KO cells (Supplementary Fig. S7E and S7F). Comparison of chromatin accessibility changes with IKZF1-bound regions in RPMI8226 indicated that 84.8% of chromatin accessibility gains directly overlapped IKZF1-bound regions (Supplementary Fig. S7G and S7H), suggesting a direct effect. Chromatin accessibility losses were primarily not at IKZF1-bound sites suggesting these changes are the result of downstream reductions in transcription factor activity. Nonetheless, the chromatin accessibility at the IGL enhancer, which is translocated in both cell types, was unaffected by lenalidomide, indicating that chromatin remains accessible (Supplementary Fig. S7I), despite the reduced potential of this key enhancer.

Given the use of a newer generation of IMiDs (also known as cereblon E3 ligase modulators or CELMoD) in MM patients, we next sought to determine if these agents could overcome the resistance associated with ETV4. Therefore, we treated RPMI8226, L363, and ARD Cas9 control and ETV4KO cells with lenalidomide, pomalidomide, and mezigdomide and found that the ETV4KO cells showed increased sensitivity to all agents when compared with Cas9 control cells (Supplementary Fig. S8A and S8B). Taken together, these data suggest that ETV4 helps overcome IKZF1 dependency in MM cells, sustains oncogenic enhancers, including those that drive MYC expression in MM, and thereby promotes resistance to IMiDs and CELMoD.

ETV4 Expression Is Associated with Poor Outcome in Primary MM Samples and Upregulated in Patients with Acquired Resistance to IMiD

We next investigated ETV4 expression in primary MM samples using RNA-seq data obtained from two myeloma trials: the MMRF CoMMpass (NCT01454297; MedRxiv 2021.08.02.21261211) and POLLUX (ref. 44; NCT02076009) studies. ETV4 expression was low to absent in most newly diagnosed (NDMM) patients but was present (≥1 FPKM) in 10.1% of the 764 NDMM samples included in the CoMMpass study (Supplementary Fig. S9A). Notably, ETV4 expression was higher in the t(11;14) Cyclin D (CD-1 and CD-2), t(4;14) MMSET (MS), and proliferation (PR) gene-expression subtypes (Supplementary Fig. S9B). Of interest, genomic copy-number alterations indicated ETV4 copy-number gains occurred in 9.3% of NDMM samples, but this was not associated with increased ETV4 expression (Supplementary Fig. S9C and S9D), suggesting that ETV4 expression primarily results from transcriptional dysregulation rather than genomic amplification. Regardless, ETV4 expression was associated with significantly reduced survival outcomes, as indicated by the PFS and overall survival (OS) Kaplan–Meier plots for CoMMpass patients treated with IMiDs as part of their first-line regimen (Fig. 7A, N = 567). Importantly, ETV4 expression remained a significant prognostic factor with different expression threshold cutoffs (Supplementary Fig. S9E) as well as when analyzed as a quantitative variable in univariate and multivariate analyses when considering other clinically relevant prognostic factors including ISS, β-2-microglobulin, age, and high-risk genetic alterations (Supplementary Table S1). Notably, when we applied a similar analysis to the patients treated with lenalidomide, and dexamethasone in the control arm of the POLLUX study, we found that ETV4 expression was also associated with a significantly worse PFS at 1 and 0.5 FPKM cutoffs (Fig. 7B; Supplementary Fig. S9F; Supplementary Data S6). In addition, we inspected ETV4 expression in the CC-4047-MM-010 (NCT01712789, hereafter MM-010) study where patients were treated with pomalidomide and dexamethasone (45). Similar to the CoMMpass and POLLUX studies, the MM-010 study showed that higher ETV4 expression corresponded with worse PFS as well as OS (Supplementary Fig. S10, Supplementary Data S6).

Figure 7.

ETV4 expression is associated with poor outcome and a proliferation gene-expression program. A, Kaplan–Meier PFS (left) and OS (right) of CoMMpass patients treated with IMiD as part of their first-line therapy (N = 567) stratified by ETV4 expression of 1 FPKM. B, PFS (left) and OS (right) analysis as in A applied to POLLUX patients treated with daratumumab, lenalidomide, and dexamethasone. C, Expression of ETV4, IKZF1, IKZF3, CRBN, IRF4, and MYC in paired samples collected at diagnosis (ND) and relapse (RR; N = 103) from CoMMpass patients (N = 47). D, Expression analysis of paired samples (N = 28) from POLLUX patients (N = 14) as in C. E, GSEA of gene sets associated with ETV4 expression in NDMM and RRMM samples from CoMMpass. Only gene sets with an FDR ≤0.05 in both analyses are shown. P values were determined using a Cox proportional hazards model Wald test (A and B) or a linear regression with a covariate for patients (C and D).

Figure 7.

ETV4 expression is associated with poor outcome and a proliferation gene-expression program. A, Kaplan–Meier PFS (left) and OS (right) of CoMMpass patients treated with IMiD as part of their first-line therapy (N = 567) stratified by ETV4 expression of 1 FPKM. B, PFS (left) and OS (right) analysis as in A applied to POLLUX patients treated with daratumumab, lenalidomide, and dexamethasone. C, Expression of ETV4, IKZF1, IKZF3, CRBN, IRF4, and MYC in paired samples collected at diagnosis (ND) and relapse (RR; N = 103) from CoMMpass patients (N = 47). D, Expression analysis of paired samples (N = 28) from POLLUX patients (N = 14) as in C. E, GSEA of gene sets associated with ETV4 expression in NDMM and RRMM samples from CoMMpass. Only gene sets with an FDR ≤0.05 in both analyses are shown. P values were determined using a Cox proportional hazards model Wald test (A and B) or a linear regression with a covariate for patients (C and D).

Close modal

Next, we investigated ETV4 expression in paired primary samples collected at diagnosis and relapse. Analysis of 76 samples obtained from 35 CoMMpass patients with paired samples revealed that ETV4 was significantly upregulated at the time of relapse (Fig. 7C), whereas no change was observed in IKZF1, IKZF3, IRF4, or MYC expression (Fig. 7C). CRBN expression was slightly downregulated (P = 0.055), which is consistent with previous reports (ref. 46; Fig. 7C). Similarly, ETV4 was found to be upregulated when we performed the same analysis on 14 patients enrolled in the POLLUX study with available paired samples (Fig. 7D). Taken together, these data suggest that ETV4 expression is associated with poor survival and upregulated in patients with resistance to IMiDs.

Finally, we sought to better understand the gene-expression program coupled with ETV4 expression in primary MM samples. We compared ETV4 expression to that of other genes in a cross-sectional analysis of 764 NDMM CoMMpass patients and identified 7,352 genes significantly associated with ETV4 expression, with roughly equal numbers of genes positively and negatively correlated with ETV4, as shown in the volcano plot (Supplementary Fig. S11A; Supplementary Data S7). Interestingly, the GSEA annotation of genes correlated with ETV4 expression indicated that the expression of cell-cycle genes, including E2F and MYC targets, was positively associated with ETV4 expression (Fig. 7E; Supplementary Data S6). Superimposing these gene sets on the volcano plot further corroborated these findings (Supplementary Fig. S11B).

Next, we determined the genes associated with ETV4 upregulation in RRMM using RNA-seq obtained from paired CoMMpass samples in conjunction with a model to control for patient-specific effects and to identify genes co-regulated with ETV4 upon relapse (see Methods). This analysis identified 28 genes that were inversely related to ETV4 expression changes and 62 genes that were positively associated with ETV4 expression changes upon relapse (Supplementary Fig. S11C and S11D). GSEA annotation of the gene-expression changes associated with ETV4 changes at relapse revealed a striking overlap with those found in the cross-sectional analysis of ETV4 expression at diagnosis, including downregulation of IL6/JAK/STAT3 signaling and inflammatory response genes and upregulation of cell cycle–related genes, including MYC and E2F targets (Fig. 7E; Supplementary Data S6). Taken together, these data implicate ETV4 expression in MM biology and demonstrate that its expression is associated with resistance to IMiDs.

In MM cells, IMiDs induce IKZF1 and IKZF3 ubiquitylation and proteasomal degradation, resulting in the transcriptional repression of MYC and IRF4 (7, 36, 37). Here, we showed that MYC transcriptional downregulation, rather than IKZF1 and IKZF3 degradation, is a better molecular correlate of IMiD efficacy. This is consistent with previous reports; however, our IKZF1 ChIP-seq data demonstrated that IMiDs via IKZF1 degradation, directly target both canonical MYC enhancers and immunoglobulin enhancers and other SEs that are juxtaposed to MYC by translocations. This is broadly relevant to MM, because MYC genetic alterations are a progression risk factor for smoldering MM (32, 34), are present in more than 35% of NDMM (29), and IMiDs are used to treat both stages of the disease (47, 48). Interestingly, we observed lenalidomide-induced downregulation of MYC only in IMiD-sensitive cells, whereas MYC expression was maintained in IMiD-resistant cells. A differentiating factor of IMiD sensitivity was the loss of chromatin-bound P300 and BRD4 coactivators from critical enhancers. Therefore, we examined other transcription factors that may maintain enhancer function and MYC expression in the absence of IKZF1 and IKZF3. This analysis implicated the ETS factor ETV4 to play a role in IMiD resistance. This was based on the similarity of the transcription factor binding motif and its expression, which correlated with IMiD resistance in both MM cell lines and primary patient samples. Using ChIP-seq, we showed that ETV4 binds to a significant portion of IKZF1-bound enhancers and that ETV4 enhancer occupancy is independent of lenalidomide-mediated IKZF1 depletion. Notably, Cas9-mediated knockout of ETV4 in three IMiD-resistant MM cell lines sensitized them to lenalidomide, pomalidomide, and mezigdomide. Together, these data demonstrate that transcriptional plasticity with the expression of extralineage TFs such as the ETS family member ETV4 sustains oncogenic enhancers in MM, overcoming IKZF1 and IKZF3 dependency, and promoting resistance to IMiD as well as CELMoD.

IMiDs have become an integral part of MM treatment (1) and have been widely used (4, 5) well before the discovery of their molecular targets (8, 9). Prior to the identification of IKZF1 and IKZF3, being major molecular targets of IMiDs, Ikzf1 and Ikzf3 were discovered to be essential for B-cell and plasma cell commitment, respectively (10, 11). However, it is still not clear how IKZF1 and/or IKZF3 regulate MYC transcription in MM cells. Our work provides new insights into the molecular basis of IKZF1 and IKZF3 functions in MM, which will help delineate the mechanisms of IMiD resistance. Nonetheless, several critical questions remain. For instance, we previously observed that MM IGL;MYC translocations are associated with reduced therapeutic benefit from IMiDs, and speculated that this may be due to the high levels of IKZF1 binding at the IGL enhancer (29). Our data suggest that it is more complex, and we cannot rule out a model whereby transcription factor combinations such as IKZF1 and ETV4 are required to maintain a threshold level of enhancer activity to maintain oncogene expression and MM survival. Additionally, previous studies (36) and those in the companion paper from Welsh and colleagues (49) found that ectopic MYC expression failed to rescue MM cells from IMiD treatment, suggesting that MYC is not the only determinant of IMiD responses. Similarly, it is possible that MYC regulates IKZF1 and IKZF3 expression, and as such MYC downregulation would further potentiate IMiD responses. Finally, we observed that IMiD treatment activated interferon response genes and downregulated MYC and cell-cycle proliferation genes and that IKZF1 and IKZF3 bound several genes in both categories. A better understanding of the cis-regulatory architecture of IKZF1 and IKZF3 may help delineate what controls the activating or repressive roles of these factors. In our analysis, we found that over half of IKZF1 sites had characteristic chromatin modifications of enhancers, consistent with the results reported in B-cell acute lymphoblastic leukemia (10). Although this alone does not provide evidence that IKZF1 functions as an activator of transcription at these sites, IMiD-mediated depletion of IKZF1 and IKZF3 resulted in the loss of P300 and BRD4 coactivators in IMiD-sensitive cells, indicating that IKZF1 is necessary for enhancer function. A striking finding of these data was IKZF1 binding to more than 75% of P300-occupied enhancers and virtually all P300-defined SEs, which are often translocated to drive oncogene expression in MM. This finding, combined with the indication that chromatin-bound P300 was lost in IMiD-responsive MM cells, provides a therapeutic rationale for the combined targeting of IKZF1/ IKZF3, and P300. This is supported by recent reports that P300 disruption effectively targets MM and other hematologic malignancies (50, 51). Indeed, in a study by Welsh and colleagues (49) IMiDs are shown to synergize with P300 inhibitors, and such rational combinations may be critical for overcoming IMiD resistance in MM.

Although IMiDs have transformed outcomes in MM, almost all patients eventually develop resistance; however, their mechanisms of resistance are not fully defined. The most common mechanisms described thus far include CRBN mutations, transcriptional downregulation, or CRBN exon 10 splicing (46, 52–54). Consistent with this, we observed a subtle reduction in CRBN expression at the time of relapse in the CoMMpass samples. These observations are further supported by CRISPR loss-of-function screens, which have identified CRBN and other cullin-RING ligase components that mediate ubiquitination as necessary for IMiD responses (55). More recently, RUNX transcription factors were shown to protect IKZF1 and IKZF3 from CRBN-mediated ubiquitination (56). However, there are few reports on how MM cells acquire resistance to IMiDs, “beyond CRBN.” As such one recent study identified upregulation of CDK6 as a targetable resistance mechanism for lenalidomide in MM (57). Here, we demonstrate that IMiD resistance is also dependent upon expression of the extralineage transcription factor ETS family member ETV4. Our genetic disruption studies indicated that ETV4 promotes IMiD resistance by sustaining some IKZF1- and IKZF3-occupied MM oncogenic enhancers, despite lenalidomide-mediated depletion of chromatin-bound IKZF1 at these loci. Of note, overexpression of ETV4 in some IMiD-sensitive MM cell lines failed to induce resistance indicating that ETV4 expression as well as a dependency on ETV4-regulated enhancers is required to induce IMiD resistance. It is also plausible that in the cell lines tested, ETV4 may not displace other chromatin-bound ETS factors or that posttranslational modifications or different binding partners may direct ETV4 away from the enhancers critical to maintaining MM cell proliferation independently of IKZF1 and IKZF3.

Of clinical relevance, we demonstrated that ETV4 expression is associated with poor survival in MM patients treated with IMiD-containing regimens from three independent trials. Whether ETV4 expression is also indicative of aggressive disease in patients treated with other regimens remains to be determined but is plausible given ETV4 binds and likely regulates immunoglobulin enhancers translocated to oncogenes in a large fraction of MM cases. Additionally, whether ETV4 regulates both intrinsic and acquired IMiD resistance is not fully understood. Given the cell lines used in this study were derived from IMiD-naïve patients as well as the observed upregulation of ETV4 in relapse samples from IMiD-treated CoMMpass and POLLUX patients, it is possible that ETV4 has a role in both forms of IMiD resistance. It will be important to better understand how often ETV4, and if other ETS, E2A, bHLH, IRF, or AP-1 factors mediate IMiD resistance. Dissecting the milieu of transcription factor occupancy and circuitry that can substitute for IKZF1 and IKZF3 will undoubtedly be a daunting task but may identify other noncanonical MM transcription factors linked to therapeutic resistance.

Cell Lines

The human MM cell lines NCI-H929, RPMI-8226, and U266 were purchased from ATCC; KMS11, MM1s, and OPM2 cells were kindly provided by Dr. Lawrence Boise (Emory University, Atlanta, GA, USA); and KMS18, JJN3, L363, XG7, XG6, and INA-6 were provided by Dr. Jonathan Keats (Translational Genomics Research Institute, Phoenix, AZ). All MM cell lines were cultured in RPMI-1640 medium (Life Technologies) containing 10% fetal bovine serum (Life Technologies), 100 U/mL penicillin, and 100 mg/mL streptomycin (Life Technologies) at 37°C under 5% CO2. INA-6, XG-7, and XG-6 were maintained in culture media supplemented with 2.5 ng/mL of human rIL6 (Sigma-Aldrich). All cell lines were authenticated by STR profiling analysis (Labcorp) and routinely tested for mycoplasma (Mycoplasma PCR Detection Kit, Sigma).

MTT Assay

Cell growth was assessed using [3-(4,5-dimethylthiazol-2-yl)-5-(3-carboxymethoxyphenyl)-2-(4-sulfophenyl)-2H-tetrazolium] according to the manufacturer's instructions (Promega).

Cell Viability Estimation by Flow Cytometry

MM cells were plated in 96-well plates at a concentration of 5 × 105 cells/mL and treated with lenalidomide (10 μmol/L), pomalidomide (1 μmol/L), and mezigdomide (0.1 μmol/L) at the indicated time points. Cells were collected, resuspended in 1× Annexin V staining buffer (BioLegend) and incubated for 10 minutes at room temperature with Annexin V–FITC antibody (BD Pharmgen) and Propidium iodide (BioVision), as per manufacturer's protocols. 300 μL PBS (Life Technologies) was added to the stained cells. Cell viability estimations were conducted using a CytoFLEX cytometer (Beckman Coulter). Based on forward and side scatter plots, 10,000 live events were recorded. The cells were displayed on in two parameter dot (density) plot (B610 on the vertical axis for PE and B525 on the horizontal axis for Annexin V–FITC) to determine the proportion of PI and Annexin V–positive cells. A cluster of cells staining strongly positive for PI and Annexin V was considered dead, and their gated percentage was used to compare target cell death among different treatment conditions. Results were shown as density plots. All flow experiment analyses were repeated three times and analyzed using Kaluza Analysis Software 2.1 (Beckman Coulter). Pairwise t tests were performed using the R function pairwise t test. A Benjamini–Hochberg (BH) adjusted P inferior to 0.05 was considered to indicate statistically significant results.

Western Blot Analysis

Equal amounts of total cellular protein from each sample were electrophoresed on sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) gels and transferred to nitrocellulose membranes. The membranes were blocked with 5% nonfat milk in Tris-buffered saline solution containing 0.1% Tween 20 for 30 minutes, probed overnight with relevant antibodies, washed, and probed with species-specific secondary antibodies coupled to horseradish peroxidase. The immunoreactive material was detected using enhanced chemiluminescence (PerkinElmer). Details of the antibodies used for the immunoblotting experiments are presented in Supplementary Table S2.

Lentiviral CRBN and ETV4 Gene Knockout

Lentiviral transduction particles (Dharmacon) were used to deliver siRNAs expressed from shRNAs to knockdown CRBN in OPM2 and KMS-11 MM cells. Lentiviral CRBN shRNA and nontarget control shRNA were produced in HEK293T packaging cells, concentrated at different MOIs, and then individually added to MM cell suspensions in the presence of 6 μg/mL polybrene and transduced for 24 hours, followed by selection with puromycin (2 μg/mL, Invitrogen) to obtain stable clones. RPMI8226, L363, and ARD cells stably expressing Cas9 (pLentiCas9-Blast, Addgene #52962) were infected with a nontargeting control sgRNAs or ETV4-targeting sgRNAs. The infected cells were then selected using 1 μg puromycin for three days. Single-cell clones were expanded and screened by immunoblotting to identify successful knockouts.

ATAC-seq

ATAC-seq was performed similar to that previously described (21, 58). Briefly, 50,000 viable cells were isolated by FACS. Nuclei were isolated by spinning at 500 g for 10 minutes, all supernatant was pipet decanted, and pelleted cells were resuspended in 25 μL Tagmentation Reaction Mix composed of 12.5 μL Tagment DNA Buffer (Illumina), 0.02% Digitonin, 0.1% Tween-20, 2.5 μL TN5 enzyme (Illumina). Tagmentation was performed for 1 hour at 37°C. Tagmented DNA was isolated by adding Tagmentation Clean-up Buffer (326 mmol/L NaCl, 109 mmol/L EDTA, 0.63% SDS, 20 μg Proteinase K) and incubating at 40°C for 1 hour and both a negative (0.7×) and positive (1.2×) SPRI bead (Kapa) selection. Libraries were amplified 12 times with Hifi polymerase (Kapa) and Nextera barcoded primers (Illumina) prior to 1× SPRI bead clean-up and sequencing on a Novaseq (Illumina).

ATAC-seq Analysis

ATAC-seq FASTQ files were quality and adapter trimmed using Trim Galore! (v0.6.4) and CutAdapt (v2.5) and mapped to the GRCh38 reference genome using bowtie2 (v2.3.5.1; ref. 59). Aligned SAM files were converted to BAM files and putative PCR duplicates were marked with Samtools (v1.10; ref. 60). Regions of chromatin accessibility were identified in each sample using MACS2 (v2.1.1.20160309; ref. 61). The union of all accessible regions was determined and reads mapping to these regions for each sample was assessed using the “summarizeOverlaps” function of the GenomicAlignments package (v1.30.0) in R (v4.1.2). Regions that overlapped ENCODE blacklisted regions were removed (62). As a quality control metric, the number of reads in peaks was determined and used to calculate a normalized accessibility score as reads per peak million (RPPM) as previously described (63) and according to the following formula:
formula

Differential chromatin accessible regions were determined using edgeR (v3.36.0; ref. 64) and imposing an FDR ≤0.01 and a fold change ≥2. Overlap of differential chromatin accessible regions with transcription factor binding motifs curated in the JASPAR database (65). Here, motifs were determined genome-wide as described in the “Transcription Factor Analysis” below, and the overlap of motifs with differentially accessible regions was determined using the Fisher exact test and FDR correction.

ChIP

For ChIP experiments, cells were harvested under the indicated conditions. ChIP experiments were performed according to standard protocols (Millipore), with minor modifications. Briefly, cells were cross-linked for 10 minutes with 1% formaldehyde at 37°C. This reaction was subsequently quenched with 125 mmol/L glycine for 5 minutes, and 5–10 million fixed cells were used for the ChIP experiments. Chromatin from fixed cells was fragmented by sonication with a Covaris E220, and the solubilized chromatin was incubated with the indicated antibody (listed in Supplementary Table S2) overnight at 4°C. Antibody–chromatin complexes were pulled down by incubation with Protein G-Dynabeads (Thermo Fisher Scientific) for 3 hours at 4°C, washed, and eluted. The samples then underwent crosslink reversal, RNase A (Roche) treatment, and proteinase K (Thermo Scientific) treatment before the captured DNA was extracted using AMP Pure beads (Beckman Coulter).

qRT-PCR Analysis

qRT-PCR was performed using TaqMan (Applied Biosystems). Expression levels were normalized to GAPDH (for gene), and the relative expression level of specific mRNA was presented as 2ΔΔCt or 2ΔCt. RNA was extracted using an RNeasy Mini Kit (Qiagen).

Molecular Data Selection

Molecular and sequencing data were obtained from the patients treated with lenalidomide at our institution. Written informed consent was obtained according to the Declaration of Helsinki and the local medical center Institutional Review Board (IRB). Samples and anonymous clinical data were obtained and analyzed using protocols approved by the local IRB.

RNA and ChIP Library Preparation and Sequencing

Sample library preparations for ChIP-seq and RNA-seq samples were performed in our laboratory using the Ion Plus Fragment Library Kit (Thermo Fisher) or Ion RNA-seq Kit (Thermo Fisher), respectively. Sequencing was performed on an Ion Proton system with a 75-bp single end for ChIP-seq and RNA-seq.

Statistical Analysis

The statistical significance of differences was determined using the Student t test. The level of significance was set at P < 0.05.

RNA-seq Analysis

RNA-seq FASTQ reads were quality and adapter trimmed using Trim Galore (v0.6.4; https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) and CutAdapt (v2.8; https://github.com/marcelm/cutadapt/), and mapped to the GRCh38 reference genome using the STAR aligner (v2.5.3a; ref. 66) using the gencode (v22) transcript database. Putative PCR duplicates were marked using Samtools (v1.13; ref. 60). Gene read counts were summarized in R (v4.1.2; ref. 67) using the “summarizeOverlaps” function in “IntersectionNotEmpty” mode of the GenomicAlignments (v1.3.0) package relative to all exons from gencode (v22) retrieved using the makeTxDbFromGFF and exonsBy functions of the GenomicFeatures package (v1.46.5; ref. 68). The reads were FPKM normalized based on the total number of autosomal and nonimmunoglobulin reads.

DEGs were determined using edgeR (v3.36.0; ref. 64) and were applied to all genes with an FPKM ≥1 in at least two samples. P values were FDR corrected using the p.adjust function, and an FDR ≤0.05 was considered significant. GSEA (v4.3.2; ref. 69) was performed using the preranked list, where rank was determined by −log10(P) × sign(fold change), and results were compared with the Hallmark gene sets (v7.5.1). Heat maps, scatter plots, and bar plots were made with bespoke code in R, available upon request.

ChIP-seq Analysis

ChIP-seq FASTQ reads were quality and adapter trimmed as described above for RNA-seq reads. ChIP-seq FASTQ files were then mapped to the GRCh38 reference genome using Bowtie2 (v2.3.5.1; ref. 59) with default parameters. Bowtie2 SAM files were converted to BAM files using Samtools (v1.13; ref. 60). ChIP-seq reads per million (RPM) normalization used total autosomal reads. ChIP-seq enrichment was determined using MACS (v2.2.7.1; ref. 70) with “-g hs -q 0.01” parameters relative to the input control for the same cell lines. ChIP-seq genome plots and heat maps were generated using the rtracklayer package (v1.54.0; ref. 71) and the bespoke R code available upon request. P300 SE analysis was performed as previously described (17), where MACS2 peaks within 15 kb were stitched together, except for those within 2.5 kb of a TSS, and the cumulative reads in these stitched regions were determined using the “countOverlaps” function of the GenomicRanges (v1.46.1) package. SEs were identified as regions with a signal above the inflection point of their enhancer rank. Differences between control and lenalidomide-treated cells were used to calculate the coverage of lenalidomide-treated samples in the regions enriched in the control. Comparison of lenalidomide-induced changes in IKZF1 binding to P300 or BRD4 restricted the analysis to regions bound by IKZF1 and P300 or BRD4 in control samples.

Transcription Factor Analysis

Transcription factor binding motif position weight matrices were obtained from the JASPAR 2020 database (65) and matched to the GRCh38 genome using the matchPWM function of the Biostrings (v2.62.0) R package with a minimum score of 0.9. Motifs enriched in RPMI8226 IKZF1 binding sites were determined relative to the same regions randomly shuffled throughout the genome 100 times. Fisher exact test was used to compute the significance P-value of motifs enriched in IKZF1-bound regions relative to shuffled regions. Only binding motifs for transcription factors expressed at ≥1 FPKM in RPMI8226 with an FDR ≤10−10 and an odds ratio ≥1.5 were considered significant. Motif logo plots were created in R, using the seqLogo (v1.60.0) package.

CoMMpass Analysis

CoMMpass structural variants, including translocations, were determined using DELLY (v0.8.7; ref. 72) as previously described (29). CoMMpass copy-number alterations were provided by TGEN, as previously described (ref. 29; MedRxiv 2021.08.02.21261211). CoMMpass RNA-seq FPKM reads were downloaded from the Genomic Data Commons. Gene-expression subtype was determined using the subtypes described by Zhan and colleagues (73). Here, a gene-expression subtype classifier for CoMMpass used RNA levels that were log2(FPKM + 1) transformed and z-score normalized prior to being multiplied by the significance analysis of the microarray score for the 50 most upregulated and 50 most downregulated genes for each subtype published by Zhan and colleagues (73). This analysis resulted in a cumulative score for each subtype for each sample, where each sample was assigned to the gene-expression subtype with the highest score.

Gene expression correlated with ETV4 expression in a cross-sectional analysis of CoMMpass NDMM samples was determined using the generalized linear model glmQLFit and glmQLFTest functions in edgeR (v3.36.0; ref. 64) according to the following formula:
formula
where GXi is the RNA-seq read count for a given gene, and ETV4 is ETV4 log2(FPKM + 1) expression value. Only genes with FPKM ≥1 in 5% of the NDMM samples were considered in the analysis. Gene expression correlated with changes in ETV4 expression in RRMM was determined using the above approach restricted to 47 patients with matched NDMM and RRMM samples and by adding a covariate for patient and a covariate for relapse as follows:
formula

Here, patient is a unique deidentified patient number, and relapse indicates whether the sample is newly diagnosed or relapsed. P values were FDR corrected, and genes with an FDR ≤0.01 were considered significant. GSEA against Hallmark gene sets was performed, as described above.

Outcome analysis of CoMMpass used interim analysis 21 (IA21) outcome data. ETV4 expression was compared with PFS and OS based on a given threshold (e.g., <1 or ≥1 FPKM) or as a quantitative variable using a Cox proportional hazards model implemented in R with the “coxph” function of the survival package (v3.4-0). Significance was determined by Wald's test with a P ≤ 0.05, which was considered significant. Multivariate analysis also considered age, stage (ISS), and high-risk genetic alterations, t(4;14), t(14;16), and del(17p), determined from whole-genome sequencing.

Data Availability Statement

Genomic data generated as part of this publication are in GEO SuperSeries GSE246436, including RNA-seq (GSE246435), ChIP-seq (GSE246429), and ATAC-seq (GSE246426). Genomic Data from the CoMMpass study (MedRxiv 2021.08.02.21261211) was obtained from the Genomic Data Commons (GDC; https://portal.gdc.cancer.gov/projects/MMRF-COMMPASS) through NIH dbGaP project phs000748.v7.p4. Access to CoMMpass was granted by the NIH dbGaP data use committee for general research use under project #14898: "Genetic and epigenetic analysis of myelomagenesis.” All CoMMpass Genomic data deposited in GDC as of October 2020 were utilized.

P. Neri reports receiving personal fees from BMS, Janssen, and Sanofi outside the submitted work. B.G. Barwick reports grants from NIH/NCI, American Society of Hematology, Paula and Rodger Riney Foundation, and Multiple Myeloma Research Foundation during the conduct of the study; personal fees from Multiple Myeloma Research and Department of Defense outside the submitted work. N. Stong reports status as employee stockholder of Bristol Myers Squibb. E. Flynt reports other support from Bristol Myers Squibb during the conduct of the study and other support from Bristol Myers Squibb outside the submitted work. J.J. Keats reports grants from Genentech and personal fees from Janssen Pharmaceuticals outside the submitted work. S. Lonial reports personal fees from BMS, GSK, Celgene, Novartis, AbbVie, Genetech, Pfizer, Regeneron, Janssen, Takeda, and Amgen during the conduct of the study; personal fees from TG Therapeutics outside the submitted work. P.L. Bergsagel reports grants from NCI during the conduct of the study; personal fees from CellCentric, Pfizer, and grants from Pfizer outside the submitted work; in addition, P.L. Bergsagel has a patent for hCRBN with royalties paid from Novartis and a patent for Vk*MYC mice with royalties paid from Pfizer. L.H. Boise reports personal fees from AbbVie and AstraZeneca outside the submitted work. N.J. Bahlis reports grants from Pfizer; personal fees from Bristol Myer Squibb, Janssen, AbbVie, Genentech, and Pfizer outside the submitted work. No disclosures were reported by the other authors.

P. Neri: Data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. B.G. Barwick: Data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. D. Jung: Investigation. J.C. Patton: Investigation. R. Maity: Investigation. I. Tagoug: Investigation. C.K. Stein: Investigation. R. Tilmont: Investigation, writing–review and editing. N. Leblay: Investigation. S. Ahn: Investigation. H. Lee: Investigation. S.J. Welsh: Investigation. D.L. Riggs: Investigation. N. Stong: Formal analysis, investigation, visualization, methodology. E. Flynt: Formal analysis, investigation. A. Thakurta: Investigation. J.J. Keats: Investigation. S. Lonial: Investigation, writing–review and editing. P.L. Bergsagel: Supervision, investigation. L.H. Boise: Formal analysis, supervision, funding acquisition, investigation, visualization, writing–original draft, writing–review and editing. N.J. Bahlis: Conceptualization, resources, formal analysis, supervision, funding acquisition, investigation, visualization, methodology, writing–review and editing.

We thank the patients enrolled in the MMRF CoMMpass and POLLUX trials and the MMRF Consortium. B.G. Barwick was supported by the ASH Scholar Award and NCI Award 1K22CA266739. L.H. Boise and B.G. Barwick were supported by the Rodger and Paula Riney Foundation and NCI award 3 R01 CA192844-04S1. B.G. Barwick, P. Neri, N.J. Bahlis, and L.H. Boise were supported by the MMRF Answer Fund Award. This work was also supported by grants from the MMRF [principal investigator (PI), N. Bahlis, co-PI: P. Neri] and Terry Fox Research Institute (M4 study; PI, N. Bahlis, co-PI, P. Neri).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Blood Cancer Discovery Online (https://bloodcancerdiscov.aacrjournals.org/).

1.
Singhal
S
,
Mehta
J
,
Desikan
R
,
Ayers
D
,
Roberson
P
,
Eddlemon
P
, et al
.
Antitumor activity of thalidomide in refractory multiple myeloma
.
N Engl J Med
1999
;
341
:
1565
71
.
2.
Raza
S
,
Safyan
RA
,
Lentzsch
S
.
Immunomodulatory drugs (IMiDs) in multiple myeloma
.
Curr Cancer Drug Targets
2017
;
17
:
846
57
.
3.
Quach
H
,
Ritchie
D
,
Stewart
AK
,
Neeson
P
,
Harrison
S
,
Smyth
MJ
, et al
.
Mechanism of action of immunomodulatory drugs (IMiDS) in multiple myeloma
.
Leukemia
2010
;
24
:
22
32
.
4.
Dimopoulos
M
,
Spencer
A
,
Attal
M
,
Prince
HM
,
Harousseau
J-L
,
Dmoszynska
A
, et al
.
Lenalidomide plus dexamethasone for relapsed or refractory multiple myeloma
.
N Engl J Med
2007
;
357
:
2123
32
.
5.
Weber
DM
,
Chen
C
,
Niesvizky
R
,
Wang
M
,
Belch
A
,
Stadtmauer
EA
, et al
.
Lenalidomide plus dexamethasone for relapsed multiple myeloma in North America
.
N Engl J Med
2007
;
357
:
2133
42
.
6.
Ito
T
,
Ando
H
,
Suzuki
T
,
Ogura
T
,
Hotta
K
,
Imamura
Y
, et al
.
Identification of a primary target of thalidomide teratogenicity
.
Science
2010
;
327
:
1345
50
.
7.
Zhu
YX
,
Braggio
E
,
Shi
C-X
,
Bruins
LA
,
Schmidt
JE
,
Van Wier
S
, et al
.
Cereblon expression is required for the antimyeloma activity of lenalidomide and pomalidomide
.
Blood
2011
;
118
:
4771
9
.
8.
Lu
G
,
Middleton
RE
,
Sun
H
,
Naniong
M
,
Ott
CJ
,
Mitsiades
CS
, et al
.
The myeloma drug lenalidomide promotes the cereblon-dependent destruction of Ikaros proteins
.
Science
2014
;
343
:
305
9
.
9.
Krönke
J
,
Udeshi
ND
,
Narla
A
,
Grauman
P
,
Hurst
SN
,
McConkey
M
, et al
.
Lenalidomide causes selective degradation of IKZF1 and IKZF3 in multiple myeloma cells
.
Science
2014
;
343
:
301
5
.
10.
Georgopoulos
K
,
Moore
DD
,
Derfler
B
.
Ikaros, an early lymphoid-specific transcription factor and a putative mediator for T cell commitment
.
Science
1992
;
258
:
808
12
.
11.
Morgan
B
,
Sun
L
,
Avitahl
N
,
Andrikopoulos
K
,
Ikeda
T
,
Gonzales
E
, et al
.
Aiolos, a lymphoid restricted transcription factor that interacts with Ikaros to regulate lymphocyte differentiation
.
EMBO J
1997
;
16
:
2004
13
.
12.
Wang
JH
,
Avitahl
N
,
Cariappa
A
,
Friedrich
C
,
Ikeda
T
,
Renold
A
, et al
.
Aiolos regulates B cell activation and maturation to effector state
.
Immunity
1998
;
9
:
543
53
.
13.
Mercola
M
,
Wang
X-F
,
Olsen
J
,
Calame
K
.
Transcriptional enhancer elements in the mouse immunoglobulin heavy chain locus
.
Science
1983
;
221
:
663
5
.
14.
Gillies
SD
,
Morrison
SL
,
Oi
VT
,
Tonegawa
S
.
A tissue-specific transcription enhancer element is located in the major intron of a rearranged immunoglobulin heavy chain gene
.
Cell
1983
;
33
:
717
28
.
15.
Banerji
J
,
Olson
L
,
Schaffner
W
.
A lymphocyte-specific cellular enhancer is located downstream of the joining region in immunoglobulin heavy chain genes
.
Cell
1983
;
33
:
729
40
.
16.
Patterson
DG
,
Kania
AK
,
Zuo
Z
,
Scharer
CD
,
Boss
JM
.
Epigenetic gene regulation in plasma cells
.
Immunol Rev
2021
;
303
:
8
22
.
17.
Hnisz
D
,
Abraham
BJ
,
Lee
TI
,
Lau
A
,
Saint-André
V
,
Sigova
AA
, et al
.
Super-enhancers in the control of cell identity and disease
.
Cell
2013
;
155
:
934
.
18.
Lovén
J
,
Hoke
HA
,
Lin
CY
,
Lau
A
,
Orlando
DA
,
Vakoc
CR
, et al
.
Selective inhibition of tumor oncogenes by disruption of super-enhancers
.
Cell
2013
;
153
:
320
34
.
19.
Jin
Y
,
Chen
K
,
De Paepe
A
,
Hellqvist
E
,
Krstic
AD
,
Metang
L
, et al
.
Active enhancer and chromatin accessibility landscapes chart the regulatory network of primary multiple myeloma
.
Blood
2018
;
131
:
2138
50
.
20.
Jia
Y
,
Zhou
J
,
Tan
TK
,
Chung
T-H
,
Wong
RWJ
,
Chooi
J-Y
, et al
.
Myeloma-specific superenhancers affect genes of biological and clinical relevance in myeloma
.
Blood Cancer J
2021
;
11
:
32
.
21.
Barwick
BG
,
Gupta
VA
,
Matulis
SM
,
Patton
JC
,
Powell
DR
,
Gu
Y
, et al
.
Chromatin accessibility identifies regulatory elements predictive of gene expression and disease outcome in multiple myeloma
.
Clin Cancer Res
2021
;
27
:
3178
89
.
22.
Hu
Y
,
Zhang
Z
,
Kashiwagi
M
,
Yoshida
T
,
Joshi
I
,
Jena
N
, et al
.
Superenhancer reprogramming drives a B-cell-epithelial transition and high-risk leukemia
.
Genes Dev
2016
;
30
:
1971
90
.
23.
Bergsagel
PL
,
Chesi
M
,
Nardini
E
,
Brents
LA
,
Kirby
SL
,
Kuehl
WM
.
Promiscuous translocations into immunoglobulin heavy chain switch regions in multiple myeloma
.
Proc Natl Acad Sci U S A
1996
;
93
:
13931
6
.
24.
Bergsagel
PL
,
Kuehl
WM
,
Zhan
F
,
Sawyer
J
,
Barlogie
B
,
Shaughnessy
J
.
Cyclin D dysregulation: an early and unifying pathogenic event in multiple myeloma
.
Blood
2005
;
106
:
296
303
.
25.
Chesi
M
,
Nardini
E
,
Brents
LA
,
Schrock
E
,
Ried
T
,
Kuehl
WM
, et al
.
Frequent translocation t(4;14)(p16.3;q32.3) in multiple myeloma is associated with increased expression and activating mutations of fibroblast growth factor receptor 3
.
Nat Genet
1997
;
16
:
260
4
.
26.
Chesi
M
,
Nardini
E
,
Lim
RSC
,
Smith
KD
,
Kuehl
WM
,
Bergsagel
PL
.
The t(4;14) translocation in myeloma dysregulates both FGFR3 and a novel gene, MMSET, resulting in IgH/MMSET hybrid transcripts
.
Blood
1998
;
92
:
3025
34
.
27.
Chesi
M
,
Bergsagel
PL
,
Shonukan
OO
,
Martelli
ML
,
Brents
LA
,
Chen
T
, et al
.
Frequent dysregulation of the c-maf proto-oncogene at 16q23 by translocation to an Ig locus in multiple myeloma
.
Blood
1998
;
91
:
4457
63
.
28.
Shou
Y
,
Martelli
ML
,
Gabrea
A
,
Qi
Y
,
Brents
LA
,
Roschke
A
, et al
.
Diverse karyotypic abnormalities of the c-myc locus associated with c-myc dysregulation and tumor progression in multiple myeloma
.
Proc Nat Acad Sci U S A
2000
;
97
:
228
33
.
29.
Barwick
BG
,
Neri
P
,
Bahlis
NJ
,
Nooka
AK
,
Dhodapkar
MV
,
Jaye
DL
, et al
.
Multiple myeloma immunoglobulin lambda translocations portend poor prognosis
.
Nat Commun
2019
;
10
:
1911
.
30.
Rustad
EH
,
Yellapantula
VD
,
Glodzik
D
,
Maclachlan
KH
,
Diamond
B
,
Boyle
EM
, et al
.
Revealing the impact of structural variants in multiple myeloma
.
Blood Cancer Discov
2020
;
1
:
258
73
.
31.
Affer
M
,
Chesi
M
,
Chen
WD
,
Keats
JJ
,
Demchenko
YN
,
Tamizhmani
K
, et al
.
Promiscuous MYC locus rearrangements hijack enhancers but mostly super-enhancers to dysregulate MYC expression in multiple myeloma
.
Leukemia
2014
;
28
:
1725
35
.
32.
Misund
K
,
Keane
N
,
Stein
CK
,
Asmann
YW
,
Day
G
,
Welsh
S
, et al
.
MYC dysregulation in the progression of multiple myeloma
.
Leukemia
2020
;
34
:
322
6
.
33.
Barwick
BG
,
Gupta
VA
,
Vertino
PM
,
Boise
LH
.
Cell of origin and genetic alterations in the pathogenesis of multiple myeloma
.
Front Immunol
2019
;
10
:
1121
.
34.
Bustoros
M
,
Sklavenitis-Pistofidis
R
,
Park
J
,
Redd
R
,
Zhitomirsky
B
,
Dunford
AJ
, et al
.
Genomic profiling of smoldering multiple myeloma identifies patients at a high risk of disease progression
.
J Clin Oncol
2020
;
38
:
2380
9
.
35.
Chng
W
,
Huang
G
,
Chung
T
,
Ng
S
,
Gonzalez-Paz
N
,
Troska-Price
T
, et al
.
Clinical and biological implications of MYC activation: a common difference between MGUS and newly diagnosed multiple myeloma
.
Leukemia
2011
;
25
:
1026
35
.
36.
Fedele
PL
,
Willis
SN
,
Liao
Y
,
Low
MS
,
Rautela
J
,
Segal
DH
, et al
.
IMiDs prime myeloma cells for daratumumab-mediated cytotoxicity through loss of Ikaros and Aiolos
.
Blood
2018
;
132
:
2166
78
.
37.
Bjorklund
CC
,
Lu
L
,
Kang
J
,
Hagner
PR
,
Havens
CG
,
Amatangelo
M
, et al
.
Rate of CRL4CRBN substrate Ikaros and Aiolos degradation underlies differential activity of lenalidomide and pomalidomide in multiple myeloma cells by regulation of c-Myc and IRF4
.
Blood Cancer J
2015
;
5
:
e354
.
38.
Sampaio
EP
,
Sarno
EN
,
Galilly
R
,
Cohn
ZA
,
Kaplan
G
.
Thalidomide selectively inhibits tumor necrosis factor alpha production by stimulated human monocytes
.
J Exp Med
1991
;
173
:
699
703
.
39.
Garcia-Diaz
A
,
Shin
DS
,
Moreno
BH
,
Saco
J
,
Escuin-Ordinas
H
,
Rodriguez
GA
, et al
.
Interferon receptor signaling pathways regulating PD-L1 and PD-L2 expression
.
Cell Rep
2017
;
19
:
1189
201
.
40.
De Matos Simoes
R
,
Shirasaki
R
,
Downey-Kopyscinski
SL
,
Matthews
GM
,
Barwick
BG
,
Gupta
VA
, et al
.
Genome-scale functional genomics identify genes preferentially essential for multiple myeloma cells compared to other neoplasias
.
Nat Cancer
2023
;
4
:
754
73
.
41.
Kaya-Okur
HS
,
Wu
SJ
,
Codomo
CA
,
Pledger
ES
,
Bryson
TD
,
Henikoff
JG
, et al
.
CUT&Tag for efficient epigenomic profiling of small samples and single cells
.
Nat Commun
2019
;
10
:
568915
.
42.
Roe
J-S
,
Mercan
F
,
Rivera
K
,
Pappin
DJ
,
Vakoc
CR
.
BET bromodomain inhibition suppresses the function of hematopoietic transcription factors in acute myeloid leukemia
.
Mol Cell
2015
;
58
:
1028
39
.
43.
Buenrostro
JD
,
Giresi
PG
,
Zaba
LC
,
Chang
HY
,
Greenleaf
WJ
.
Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position
.
Nat Methods
2013
;
10
:
1213
8
.
44.
Dimopoulos
MA
,
Oriol
A
,
Nahi
H
,
San-Miguel
J
,
Bahlis
NJ
,
Usmani
SZ
, et al
.
Daratumumab, lenalidomide, and dexamethasone for multiple myeloma
.
N Engl J Med
2016
;
375
:
1319
31
.
45.
Dimopoulos
MA
,
Palumbo
A
,
Corradini
P
,
Cavo
M
,
Delforge
M
,
Di Raimondo
F
, et al
.
Safety and efficacy of pomalidomide plus low-dose dexamethasone in STRATUS (MM-010): a phase 3b study in refractory multiple myeloma
.
Blood
2016
;
128
:
497
503
.
46.
Gooding
S
,
Ansari-Pour
N
,
Towfic
F
,
Estévez
MO
,
Chamberlain
PP
,
Tsai
KT
, et al
.
Multiple cereblon genetic changes are associated with acquired resistance to lenalidomide or pomalidomide in multiple myeloma
.
Blood
2021
;
137
:
232
7
.
47.
Mateos
M-V
,
Hernández
M-T
,
Giraldo
P
,
de la Rubia
J
,
de Arriba
F
,
Corral
LL
, et al
.
Lenalidomide plus dexamethasone for high-risk smoldering multiple myeloma
.
N Engl J Med
2013
;
369
:
438
47
.
48.
Lonial
S
,
Jacobus
S
,
Fonseca
R
,
Weiss
M
,
Kumar
S
,
Orlowski
RZ
, et al
.
Randomized trial of lenalidomide versus observation in smoldering multiple myeloma
.
J Clin Oncol;
2020
;
38
:
1126
37
.
49.
Welsh
SJ
,
Barwick
BG
,
Meermeier
EW
,
Riggs
DL
,
Shi
C-X
,
Zhu
YX
, et al
.
Transcriptional heterogeneity overcomes super-enhancer disrupting drug combinations in multiple myeloma
.
Blood Cancer Discov
2024
;
5
:
34
55
.
50.
Hogg
SJ
,
Motorna
O
,
Cluse
LA
,
Johanson
TM
,
Coughlan
HD
,
Raviram
R
, et al
.
Targeting histone acetylation dynamics and oncogenic transcription by catalytic P300/CBP inhibition
.
Mol Cell
2021
;
81
:
2183
200
.
51.
Vannam
R
,
Sayilgan
J
,
Ojeda
S
,
Karakyriakou
B
,
Hu
E
,
Kreuzer
J
, et al
.
Targeted degradation of the enhancer lysine acetyltransferases CBP and p300
.
Cell Chem Biol
2021
;
28
:
503
14
.
52.
Kortüm
KM
,
Mai
EK
,
Hanafiah
NH
,
Shi
C-X
,
Zhu
Y-X
,
Bruins
L
, et al
.
Targeted sequencing of refractory myeloma reveals a high incidence of mutations in CRBN and Ras pathway genes
.
Blood
2016
;
128
:
1226
33
.
53.
Huang
S-Y
,
Lin
C-W
,
Lin
H-H
,
Yao
M
,
Tang
J-L
,
Wu
S-J
, et al
.
Expression of cereblon protein assessed by immunohistochemicalstaining in myeloma cells is associated with superior response of thalidomide- and lenalidomide-based treatment, but not bortezomib-based treatment, in patients with multiple myeloma
.
Ann Hematol
2014
;
93
:
1371
80
.
54.
Broyl
A
,
Kuiper
R
,
van Duin
M
,
van der Holt
B
,
el Jarari
L
,
Bertsch
U
, et al
.
High cereblon expression is associated with better survival in patients with newly diagnosed multiple myeloma treated with thalidomide maintenance
.
Blood
2013
;
121
:
624
7
.
55.
Sievers
QL
,
Gasser
JA
,
Cowley
GS
,
Fischer
ES
,
Ebert
BL
.
Genome-wide screen identifies cullin-RING ligase machinery required for lenalidomide-dependent CRL4CRBN activity
.
Blood
2018
;
132
:
1293
303
.
56.
Zhou
N
,
Gutierrez-Uzquiza
A
,
Zheng
XY
,
Chang
R
,
Vogl
DT
,
Garfall
AL
, et al
.
RUNX proteins desensitize multiple myeloma to lenalidomide via protecting IKZFs from degradation
.
Leukemia
2019
;
33
:
2006
21
.
57.
Ng
YLD
,
Ramberger
E
,
Bohl
SR
,
Dolnik
A
,
Steinebach
C
,
Conrad
T
, et al
.
Proteomic profiling reveals CDK6 upregulation as a targetable resistance mechanism for lenalidomide in multiple myeloma
.
Nat Commun
2022
;
13
:
1009
.
58.
Corces
MR
,
Buenrostro
JD
,
Wu
B
,
Greenside
PG
,
Chan
SM
,
Koenig
JL
, et al
.
Lineage-specific and single-cell chromatin accessibility charts human hematopoiesis and leukemia evolution
.
Nat Genet
2016
;
48
:
1193
203
.
59.
Langmead
B
,
Salzberg
SL
.
Fast gapped-read alignment with Bowtie 2
.
Nat Methods
2012
;
9
:
357
9
.
60.
Li
H
,
Handsaker
B
,
Wysoker
A
,
Fennell
T
,
Ruan
J
,
Homer
N
, et al
.
The sequence alignment/map format and SAMtools
.
Bioinformatics
2009
;
25
:
2078
9
.
61.
Feng
J
,
Liu
T
,
Qin
B
,
Zhang
Y
,
Liu
XS
.
Identifying ChIP-seq enrichment using MACS
.
Nat Protoc
2012
;
7
:
1728
40
.
62.
Amemiya
HM
,
Kundaje
A
,
Boyle
AP
.
The ENCODE blacklist: Identification of problematic regions of the genome
.
Sci Rep
2019
;
9
:
1
5
.
63.
Barwick
BG
,
Scharer
CD
,
Martinez
RJ
,
Price
MJ
,
Wein
AN
,
Haines
RR
, et al
.
B cell activation and plasma cell differentiation are inhibited by de novo DNA methylation
.
Nat Commun
2018
;
9
:
1900
.
64.
Robinson
MD
,
McCarthy
DJ
,
Smyth
GK
.
edgeR: a Bioconductor package for differential expression analysis of digital gene expression data
.
Bioinformatics
2010
;
26
:
139
40
.
65.
Sandelin
A
.
JASPAR: an open-access database for eukaryotic transcription factor binding profiles
.
Nucleic Acids Res
2004
;
32
:
91D
94
.
66.
Dobin
A
,
Davis
CA
,
Schlesinger
F
,
Drenkow
J
,
Zaleski
C
,
Jha
S
, et al
.
STAR: ultrafast universal RNA-seq aligner
.
Bioinformatics
2013
;
29
:
15
21
.
67.
Gentleman
RC
,
Carey
VJ
,
Bates
DM
,
Bolstad
B
,
Dettling
M
,
Dudoit
S
, et al
.
Bioconductor: open software development for computational biology and bioinformatics
.
Genome Biol
2004
;
5
:
R80
.
68.
Lawrence
M
,
Huber
W
,
Pagès
H
,
Aboyoun
P
,
Carlson
M
,
Gentleman
R
, et al
.
Software for computing and annotating genomic ranges
.
PLoS Comput Biol
2013
;
9
:
e1003118
.
69.
Subramanian
A
,
Tamayo
P
,
Mootha
VK
,
Mukherjee
S
,
Ebert
BL
,
Gillette
MA
, et al
.
Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles
.
Proc Nat Acad Sci U S A
2005
;
102
:
15545
50
.
70.
Zhang
Y
,
Liu
T
,
Meyer
CA
,
Eeckhoute
J
,
Johnson
DS
,
Bernstein
BE
, et al
.
Model-based analysis of ChIP-seq (MACS)
.
Genome Biol
2008
;
9
:
R137
.
71.
Lawrence
M
,
Gentleman
R
,
Carey
V
.
rtracklayer: an R package for interfacing with genome browsers
.
Bioinformatics
2009
;
25
:
1841
2
.
72.
Rausch
T
,
Zichner
T
,
Schlattl
A
,
Stütz
AM
,
Benes
V
,
Korbel
JO
.
DELLY: structural variant discovery by integrated paired-end and split-read analysis
.
Bioinformatics
2012
;
28
:
i333
9
.
73.
Zhan
F
,
Huang
Y
,
Colla
S
,
Stewart
JP
,
Hanamura
I
,
Gupta
S
, et al
.
The molecular classification of multiple myeloma
.
Blood
2006
;
108
:
2020
8
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.