A multiomics approach reveals RNA dynamics promote cellular sensitivity to DNA hypomethylation

The search for new approaches in cancer therapy requires a mechanistic understanding of cancer vulnerabilities and anti-cancer drug mechanisms of action. Problematically, some effective therapeutics target cancer vulnerabilities that have poorly defined mechanisms of anti-cancer activity. One such drug is decitabine, a frontline therapeutic approved for the treatment of high-risk acute myeloid leukemia (AML). Decitabine is thought to kill cancer cells selectively via inhibition of DNA methyltransferase enzymes, but the genes and mechanisms involved remain unclear. Here, we apply an integrated multiomics and CRISPR functional genomics approach to identify genes and processes associated with response to decitabine in AML cells. Our integrated multiomics approach reveals RNA dynamics are key regulators of DNA hypomethylation induced cell death. Specifically, regulation of RNA decapping, splicing and RNA methylation emerge as critical regulators of cellular response to decitabine.

role in determining clinical outcomes following treatment with decitabine. This result suggests that 1 decitabine's anti-cancer activity in AML occurs through a TP53 independent mechanism. Given the 2 central role TP53 plays in canonical apoptotic BCL2 family protein dependent programmed cell 3 death, at one level this study appears to contradict recent clinical trial results in AML which 4 demonstrated superior clinical outcomes from the combination of HMA and venetoclax, a BCL2 5 inhibitor thought to drive programmed cell death in cancer cells 9 . One explanation that could explain 6 both sets of clinical observations is that HMAs may drive cell death via an unknown TP53 7 independent apoptotic pathway. A more robust understanding of decitabine's mechanisms of anti-8 cancer activity in TP53-mutant tumors could enable innovative therapeutic strategies and a better 9 understanding of patients who do and do not respond robustly to HMAs. An alternate hypothesis for 10 how HMAs kill cancer cells arises from the observation that treatment with HMAs results in 11 accumulation of non-canonical transcripts including inverted SINE elements, endogenous retroviral 12 elements and cryptic transcription start sites encoded in long terminal repeats which collectively act 13 to induce immune activation 10-14 . Lastly, it has also been suggested that HMAs induce cellular 14 differentiation in AML which may contribute to therapeutic efficacy 15 . 15 16 To identify genes that modulate decitabine's anti-cancer activity in high-risk AML in an unbiased 17 manner, we performed genome-scale CRISPR genetic screens and integrated this data with 18 multiomics measurements of decitabine response. Our results recapitulate multiple known factors 19 which modulate response to decitabine, including DCK, SLC29A1, MCL1 and BCL2, indicating the 20 utility and robustness of our approach for interrogating the biology of decitabine in AML 9, [16][17][18][19][20][21][22] . 21 Central to our study was the finding that epitranscriptomic RNA modification and RNA quality 22 control pathways effectively modulate response to decitabine in AML. In short, biologically we have 23 identified regulatory connections between DNA methylation, RNA methylation and RNA quality 24 control pathways, and translationally we have nominated unexpected and clinically actionable 25 biomarkers and combination therapies which in the future may be used to potentiate decitabine's 26 anti-cancer activity. 27 28 time-zero samples, we also collected samples after growing the library in the presence and absence 1 of decitabine (in biological duplicates). Next-generation sequencing was used to quantify the relative 2 abundance of cells expressing each sgRNA in each sample. We then used measurements across the 3 entire library to calculate sgRNA-and gene-level phenotypic scores ( Figure S2A). For each cell line, 4 results obtained from the replicate screens were highly correlated with high data quality in both the 5 DMSO and decitabine experiments (Figures S2B-E; Table S1). Analysis of our decitabine screen 6 data revealed a large number of genes that modulate cellular response to decitabine (1293 genes with 7 Mann-Whitney p-value < 0.05 and absolute value of rho score > 0.1) ( Figure 1C; Table S1). These 8 results may reflect the pleiotropic nature of DNA methylation biology. 9 10 Diverse biological processes including RNA dynamics and apoptosis modulate AML cellular 11 response to decitabine 12 13 Initial inspection of top hits from our decitabine CRISPRi screen in HL-60 cells recapitulated a 14 number of genes whose knockdown is known to impact drug resistance and sensitivity ( Figure 1C; 15 Table S1). For example, the top resistance hit was DCK, which phosphorylates decitabine resulting 16 in conversion of the pro-drug to the active drug. Another top resistance hit was SLC29A1, which is a 17 solute carrier protein required for decitabine entry into cells. Lastly, DCTD is thought to play a role 18 in the metabolism of decitabine and is revealed as a strong resistance hit as well 29 . We also observed 19 that knockdown of BCL2 and MCL1 sensitizes HL-60i to decitabine, as expected from the clinical 20 literature which suggests decitabine induces BCL2 family protein mediated cell death 20,21,30 . The 21 recapitulation of known positive control hits in our screens indicate the utility and robustness of our 22 approach for interrogating the biology of decitabine in AML.

24
Buoyed by these positive endogenous controls, we next examined the remaining CRISPRi hits to 25 search for new biological insights and to generate hypotheses on the cellular mechanisms of 26 decitabine action. First, we noted that the pathway-level analysis of our screen identifies mRNA 27 processing pathways as a top-scoring enriched term ( Figure S2F; Table S2) [31][32][33] . Further manual 28 analysis of these top hits revealed a strong enrichment for two specific RNA biological processes 29 (Tables S1-S2). Specifically, we observed that repression of RNA decapping enzymes such as 30 DCP1A, DCP2 and DCPS sensitizes HL-60 to decitabine ( Figure 1C; Table S1). We also observed 31 that repression of multiple genes (METTL3, YTHDF2, ZC3H13 and CBLL1) that regulate RNA 32 methylation marks, specifically N 6 -methyladenosine or m 6 A promoted resistance to decitabine 33 ( Figure 1C; Table S1). Together, these observations suggest that modulation of specific RNA 34 regulatory pathways is a key determinant of response to DNA hypomethylation induced by 35 decitabine. To independently validate the results from our screen, we chose 10 hit genes from our 36 decitabine HL-60 CRISPRi screen (2 sgRNAs/gene) and used a mixed competition fluorescence cell 37 survival CRISPRi knockdown assay to measure how perturbation of individual genes modulates 38 response to decitabine. Our validation experiments demonstrated the reproducibility of our CRISRPi 39 genome-scale screen measurements across all the resistance and sensitivity genes tested ( Figures 1D-40 F). Interestingly, we observed that repression of PTEN, a tumor suppressor gene that is commonly 41 mutated in cancer, sensitized HL-60 cells to decitabine ( Figure 1E). 42 1 Comparative CRISPRi functional genomics experiments reveal common and specific genes 2 modulating cellular response to decitabine in additional AML models 3 4 Given the known heterogeneity of AML, we chose to perform genome-scale CRISPRi screens in two 5 additional AML models to further examine the degree of common and specific mechanisms across 6 cell lines that regulate cellular response to decitabine. For this we used SKM-1 and MOLM-13 cells, 7 which are established models of AML. Comparing the known driver mutations in these AML models 8 we noted that SKM-1 is TP53 and KRAS mutant, which similarly to HL-60, captures the biology of 9 high-risk AML and more generally of an aggressive human cancer. Meanwhile, MOLM-13 is FLT3-10 ITD and MLL-fusion but TP53 wild-type. 11 12 We engineered CRISPRi cell lines for each model and performed genome-scale CRISPRi screens to 13 identify genes that regulate response to decitabine (~IC30; 15-100 nM) as described above and 14 compared the results with the HL-60 screen ( Figure S3; Table S3). Similar to the HL-60 screen, we 15 observed that the SKM-1 and MOLM-13screens also captured mRNA processing as an enriched 16 term across top hits(Figures S3C and S3F; Table S4) and positive control genes whose knockdown is 17 known to impact drug resistance, namely DCK, SLC29A1 and DCTD (Figures S3G-3H; Table S3). 18 As expected from the heterogeneity of AML, we also observed differences across cell lines with 19 respect to genes that modulate response to decitabine (Figures S3G-H; Table S3). Interestingly, the 20 two cell lines classified as TP53-inactive (HL-60 and SKM-1), and are representative of the high-risk 21 AML patient cohort that benefits from the combination therapy of decitabine and venetoclax, 22 revealed BCL2 and MCL1 as sensitizing hits in the presence of decitabine, while the TP53-wild-type 23 cell line (MOLM-13) did not.

25
Additionally, we observed that repression of METTL3, a gene that regulates RNA methylation 26 marks, specifically N6-methyladenosine or m 6 A, promoted resistance to decitabine across all three 27 cell lines (Figures S3G-H; Table S3). In contrast, repression of genes encoding RNA decapping 28 enzymes such as DCP2 and DCPS sensitized HL-60 and SKM-1 cells, but not MOLM-13 cells, to 29 decitabine treatment. In summary, comparison of genome-scale decitabine CRISPRi screens in 3 30 AML models reveals common and unique regulators of response. This finding is in line with our 31 understanding of the heterogeneity of AML biology and suggests that therapeutic strategies in AML 32 should be evaluated in multiple models representative of diverse tumors. 33 34 RNA decapping modulates response to DNA hypomethylation induced by decitabine in AML 35 cells 36 37 We were intrigued by the connection between decitabine and RNA decapping quality control 38 processes. To begin, we confirmed that repression of DCP2 sensitizes cells to decitabine ( Figure 1E). 39 We further tested that RNA decapping is a pro-survival dependency by combining RG3039, a 40 clinical grade chemical inhibitor of DCPS, with decitabine 34-36 . We observed the combination of 41 decitabine and RG3039 had synergistic anti-cancer activity in two independent AML cell models 42 ( Figures S4A-B). We also observed that the combination of decitabine and RG3039 synergistically 1 induced caspase 3/7 activation and cell cycle arrest in HL-60 ( Figures 1G-H). Lastly, we profiled the 2 transcriptional consequences of treating cells with DMSO, decitabine alone, RG3039 alone or 3 decitabine and RG3039 together. Because previous literature has demonstrated HMAs can induce 4 expression of endogenous retroviral elements, we mapped both protein coding transcript expression 5 and ERV transcript expression. We observed that treatment with decitabine or RG3039 alone drives 6 a transcriptional response ( Figures S4C-D), and that the combination of decitabine with RG3039 7 induces transcriptional responses shared with the single drug conditions but also drug combination 8 specific transcriptional changes ( Figures S4C-D). Gene ontology analysis comparing decitabine to 9 decitabine plus RG3039 or DMSO to decitabine plus RG3039 demonstrated up regulation of term 10 enrichment for biological processes such as myeloid differentiation and immune function, as well as 11 down regulation for biological processes relating to methylation and protein translation ( Figure S4E). 12 For example, we observed the upregulation of positive regulators of TNFα cytokine production 13 specifically in the decitabine plus RG3039 condition relative to decitabine alone ( Figure 1I). Prior 14 studies have shown decitabine treatment alone can induce expression of atypical transcripts which in 15 turn can induce an inflammatory response 10,36 . Our analysis of ERV transcriptional changes 16 demonstrated that the combination of decitabine plus RG3039 strongly induced specific ERV 17 transcripts, such as LTR67B (chr6:36350628−36351191), relative to DMSO or each single drug 18 alone ( Figure 1J). Notably, most ERVs do not change expression, and changes in expression are 19 often not concordant across families or classes of ERVs ( Figures S4F-G). Together, this data 20 suggests that a possible combination therapy consisting either of two drugs (decitabine and an RNA 21 decapping inhibitor) or three drugs (decitabine, venetoclax and a RNA decapping inhibitor) could 22 have synergistic activity for the treatment of AML.

24
Decitabine treatment of AML cells results in m 6 A hypermethylation of mRNA transcripts 25 26 As highlighted above, we observed that repression of multiple genes encoding m 6 A methylation 27 machinery promotes cellular resistance to decitabine (Figures 1C, 1F and S3G-H). Top hits included 28 the m 6 A-writer METTL3, the m 6 A-reader YTHDF2 and the methyltransferase complex components 29 ZC3H13 and CBLL1. We validated that repression of METTL3, YTHDF2, and ZC3H13 promotes 30 resistance to DNMT inhibition by decitabine treatment in HL-60i over a time course using a mixed 31 competition fluorescence cell survival CRISPRi knockdown assay ( Figure 2A). This result suggests 32 regulation of RNA methylation modulates AML cell survival upon treatment with decitabine. 33 34 To systematically examine the molecular effect of decitabine treatment on m 6 A RNA methylation, 35 we next performed methylated RNA immunoprecipitation sequencing (MeRIP-seq), a method for 36 detection of m 6 A modifications ( Figure S5A) 37 . To assess the quality of this dataset, we first 37 performed peak calling in control DMSO-treated samples followed by downstream analysis to 38 recapitulate known features of the RNA modification sites across the transcriptome. We also 39 confirmed the preferential localization of RNA methylation peaks near the stop codon, which is 40 consistent with prior literature ( Figure 2B) 38 . Lastly, we performed a motif-enrichment analysis to 41 ensure the enrichment of the RGAC ([AG]GAC) motif sequence, a known m 6 A motif, among 1 predicted peaks ( Figure 2C) 39-40 . 2 3 To then identify decitabine-induced hyper-and hypomethylated sites, we performed differential 4 RNA methylation analysis to compare treatment with decitabine to DMSO controls 41 . Interestingly, 5 we observed a significant increase in m 6 A RNA methylation peaks across mRNAs of protein coding 6 genes upon decitabine treatment ( Figure 2D; Table S5). Specifically, our analysis identified 2064 7 decitabine induced hypermethylated peaks (logFC >1 and p-value <0.005) but only 1399 8 hypomethylated peaks (logFC <−1 and p-value <0.005) ( Figures S5B-C). 9 10 It has been observed in AML cell lines and patient data that treatment with different HMAs such as 11 decitabine induces transcriptional upregulation of different ERVs including retroposons, LINEs and 12 SINEs 12,42,43 . It has also been shown that m 6 A RNA methylation regulates the levels of ERVs 44 . To 13 evaluate the effect of decitabine treatment on ERV RNA methylation, we mapped our MeRIP-seq 14 data to relevant annotations and followed similar analyses as discussed above to examine differential 15 RNA methylation changes in ERVs 45,46 . Interestingly we observed a significant enrichment of m 6 A 16 methylation peaks across retroposon, LINE and SINE transcripts upon decitabine treatment ( Figure  17 2E-F). Specifically, our analysis here identified 37, 180 and 131 hypermethylated peaks (logFC >1 18 and p-value <0.005) but only 9, 45 and 48 hypomethylated peaks (logFC <−1 and p-value <0.005) 19 for retroposon, LINE and SINE transcripts, respectively. 20 21 Taken together, our findings suggest that treatment of AML cells with decitabine results in global 22 CpG DNA hypomethylation along with a concomitant increase in m 6 A RNA methylation, and that 23 HMA anti-cancer activity in AML cells may be modulated by genes that regulate m 6 A RNA 24 methylation.

26
A multiomics approach identifies genes regulated through m 6 A modifications 27 28 RNA methylation has been implicated in various aspects of the RNA life cycle in the cell, from RNA 29 processing to RNA stability to translation, and more recently, crosstalk between epitranscriptome 30 and epigenome [47][48][49][50][51][52][53][54][55] . To further understand the connection between global DNA hypomethylation and 31 RNA dynamics in AML cells, we set out to interrogate, via an integrated multiomics approach, the 32 effects of decitabine-induced RNA hypermethylation on AML cells. Here, we aimed to integrate 33 comparisons between treatment with decitabine or DMSO from the following datasets: RNA-seq for 34 differential gene expression and RNA stability, MeRIP-seq for RNA methylation, Ribo-seq for 35 protein translation efficiency, and genome-scale CRISPRi functional genomics screening data. We 36 first performed an RNA-seq time course experiment in the HL-60 AML model ( Figure S6A) at 6, 72 37 and 120 hours following treatment with decitabine or DMSO. We used this data to perform 38 differential gene expression analysis across conditions. We also used REMBRANDTS, a method we 39 have previously developed for differential RNA stability analysis, to estimate post-transcriptional 40 modulations in relative RNA decay rates ( Figure 3A-B) 56-61 . We performed gene set enrichment 41 analysis of differential mRNA stability and expression across all three time points for the HL-60 cell 42 line ( Figures S6B-C) 62 . For expression, we observed enrichment for largely expected ontologies, 1 such as immune receptor activity and regulation of cell killing 10, 12,14,20,36 . Interestingly, for post-2 transcriptional modulations in RNA stability, we observed previously unexplored terms, such as 3 sterol biology. Moreover, to also capture patient heterogeneity, we performed RNA-seq on a panel of 4 five additional AML cell lines treated with decitabine or DMSO. Across all six AML cell lines, we 5 observed that decitabine treatment induced widespread changes in RNA transcript abundance and 6 RNA stability with varying degrees of concordant RNA expression and stability changes ( Given that RNA m 6 A methylation marks have been previously implicated in translational control, we 10 used Ribo-seq to measure changes in the translational efficiency landscape of HL-60 cells treated 11 with decitabine or DMSO 50,63 . Treatment with decitabine had little effect on translation efficiency, 12 and we did not observe a concerted change in the translation efficiency of hypermethylated mRNAs 13 ( Figure S7). In other words, changes in translation efficiency of mRNAs that are differentially 14 methylated in decitabine-treated cells are not likely to be responsible for cellular sensitivity to this 15 drug.

17
Having ruled out translational control as the mechanism through which RNA methylation may be 18 involved, we next sought to identify genes whose RNA hypermethylation drives cellular sensitivity 19 to decitabine through other post-transcriptional regulatory programs. Since m 6 A RNA methylation is 20 shown to reduce RNA stability and expression, we intersected our set of decitabine-induced 21 hypermethylated genes with those that are downregulated in decitabine treated cells, and their lower 22 expression is associated with higher sensitivity to decitabine in our functional CRISPRi screen 23 ( Figure 4A) 64,65 . In this analysis, we observed ten genes that were sensitizing hits in the CRISPRi 24 screen and upon decitabine treatment, showed RNA hypermethylated peaks and lower mRNA levels 25 ( Figure 4B). We observed that these genes collectively regulate nuclear processes (INTS5, INO80D,  26 ZNF777, MYBBP1A, RNF126, RBM14-RBM4)) or metabolism (SQLE, DHODH, PMPCA, SLC7A6). 27 To extend our observations, we also identified genes that (i) were downregulated upon decitabine 28 treatment across our panel of six AML cell lines, (ii) sensitizing hits in our HL-60 CRISPRi screen, 29 and (iii) showed hypermethylated peaks upon decitabine treatment in our MeRIP-seq HL-60 data 30 ( Figures 4C-D). Although this analysis converges on a very small number of genes, we were 31 nevertheless intrigued by the possibility that several nominated genes could serve as a link between 32 RNA methylation and the cell death induced by decitabine.

34
From this list we selected SQLE and INTS5 and validated that their mRNA abundance is decreased 35 and m 6 A methylation is increased in HL-60 cells following treatment with decitabine (Figures 4E  36 and S8A-B). Consistently, we observed that SQLE and INTS5 pre-mRNA levels do not change, 37 showing that the decreased mRNA levels are not due to a decrease in transcription. These results 38 suggest we have identified a small number of mRNAs that are likely downregulated due to increased 39 m 6 A methylation following decitabine treatment and that these same genes are functionally 40 important for cellular response to decitabine. 41 42

2
Our experiments identify previously known and unknown genes and pathways that modulate cellular 3 response to decitabine, a clinically approved HMA with poorly understood cellular mechanisms of 4 action. Our results unexpectedly reveal a key role for RNA dynamics in modulating the response to 5 DNA hypomethylation induced by decitabine. 6 7 Specifically, we observed that genes which are thought to regulate mRNA decapping promote 8 cellular resistance to decitabine. One hypothesis for why loss of RNA decapping enzyme activity 9 sensitizes AML cells to decitabine is that this RNA quality control pathway becomes an induced 10 dependency upon decitabine treatment due to repressed or aberrant transcripts that accumulate upon 11 decitabine-induced DNA hypomethylation. Alternatively, some RNA decapping proteins are also 12 key regulators of splicing, so it may be that this biology is more complex with respect to 13 transcription than currently appreciated 66,67 . 14 15 We also observed that genes which are thought to both write and read m 6 A RNA methyl marks 16 promote sensitivity to decitabine. Although many functions of m 6 A have been identified, here we 17 demonstrate a drug phenotype associated with m 6 A biology. Given all known human 18 methyltransferase enzymes use S-adenosyl methionine (SAM) as a cofactor for transfer of methyl 19 groups, an intriguing hypothesis arises in which treatment of cells with decitabine results in global 20 inhibition of DNMTs, resulting in increased SAM levels and subsequently hypermethylation of 21 mRNAs leading to cell death. To our knowledge, cross talk between methyl transferase enzymes 22 with different macromolecular substrates is not known, and this hypothesis may merit further 23 investigation.

25
Lastly, our efforts may have several translational implications for AML patients who are treated with 26 decitabine. We experimentally confirm that decitabine induces TP53 independent apoptosis in 27 experimental models. In line with this, our results genetically re-nominate a clinically efficacious 28 combination therapy of decitabine and a BCL2 inhibitor, which together likely induces synergistic 29 apoptosis. We also demonstrate in AML cellular models that the combination of decitabine and 30 RG3039, an RNA decapping inhibitor, could be a potential new route for the potentiation of HMAs 31 in AML treatment.

33
We anticipate that our study serves as an integrated multiomics resource for understanding AML 34 cellular response to decitabine and we hope these results also provide broader insight into HMAs as 35 anti-cancer drugs and aid in the design of future therapeutic strategies.

37
Acknowledgments 38 39 We would like to thank members of the Goodarzi and Gilbert labs for helpful discussions and 40 assistance, including Dr. Gary Wang, Dr. Becky Xu Hua Fu, Dr. Benedict Choi, Laine Goudy, and 41 Tanvi Joshi. We also thank Dr. Felix Krueger for helping A.A. troubleshoot Bismark. L.A.G. is 42 cells on a BD FACS Aria III. Cells were then assayed for CRISPRi activity by transducing sgRNAs 1 targeting two essential genes (PLK1, HSPA9) and assessing for relative depletion of GFP (i.e., 2 sgRNA positive cells) via flow cytometry between day 3 and day 9 post-transfection.

4
CRISPRi screen experimental procedure 5 Genome-scale CRISPRi screens were performed similarly to those previously described 23 . The 6 human CRISPRi-v2 sgRNA library (top 5 sgRNAs per gene) was transduced into HL-60, SKM-1 7 and MOLM-13 cells at 250-500-fold coverage. Cells were resuspended in lentiviral supernatant with 8 8 μg/mL polybrene in 6-well plates and centrifuged at 1000 g for 2 hours at room temperature. Cells 9 were resuspended into fresh medium following spinfection. 72 hours following infection, cells were 10 seeded at 1,000,000 cells/mL for puromycin selection (0.5-1 ug/mL). Following puromycin 11 selection, "time-zero" samples were harvested at 500x library coverage. CRISPRi screen computational analysis 24 Sequencing reads were trimmed, aligned to the human CRISPRi-v2 sgRNA library and counted 25 using a previously described pipeline (https://github.com/mhorlbeck/ScreenProcessing). Growth (γ) 26 and drug sensitivity/resistance (ρ) phenotypes were calculated based on sgRNA frequencies across 27 conditions 23,70 . Gene phenotypes were calculated by taking the mean of the top three sgRNA 28 phenotypes per gene by magnitude. Gene phenotype p-values were calculated using the Mann-29 Whitney test comparing the gene-targeting sgRNAs with a set of non-targeting control sgRNAs. For 30 genes with multiple annotated transcription start sites (TSS), sgRNAs were first clustered by TSS, 31 and the TSS with the smallest Mann-Whitney p-value was used to represent the gene. Hits were 32 defined as genes with a phenotype Z-score greater or equal to 6. Z-scores were calculated by 33 dividing the gene phenotype by the standard deviation of the non-targeting sgRNA phenotypes 23 .
(2) To identify broader pathways associated with drug 3 response irrespective of ρ phenotype direction, we performed GSEA analysis on genes ranked by 1 -4 Mann-Whitney p-value (calculated for each ρ phenotype as above) and set a minimum threshold for 5 gene set size (i.e., `min_size=200`). Then we only report positive normalized enrichment scores 6 (NES) that represent gene sets enriched by more significant Mann-Whitney p-value for the gene set's 7 ρ phenotype. 8 9 Individual sgRNA validation 10 11 Individual sgRNAs were validated using a competitive growth assay as previously described 23 . 12 Briefly, sgRNA protospacers with flanking BstXI and BlpI restriction sites were cloned into the 13 BstXI/BlpI-digested pCRISPRia-v2 plasmid (Addgene #84832). Protospacer sequences are listed in 14 Cell viability assay and Bliss excess score calculation 2 Cells were seeded into 96-well plates at 100,000 cells/mL in duplicate and were treated with 3 decitabine (seven-point 1:3 dilution series from 0.5 uM to 0.002 uM), RG3039 (seven-point 1:4 4 dilution series from 10 uM to 0.010 uM) or the combination of both drugs at all possible dose 5 combinations. Control cells treated with DMSO were counted at day 3, and all cells were split at the 6 ratio required to dilute control cells to a concentration of 100,000 cells/mL. Raw fluorescence units 7 (RLUs) were assessed at day 3, day 5 and day 7 for each condition using the CellTiter-Glo (CTG) 8 luminescence-based assay (Promega). Diluted CTG reagent (100 uL 1:4 CTG reagent to PBS) was 9 added to cells (100 uL) and the mixture was pipetted up and down to ensure complete cell lysis.

10
Luminescence was then assayed using a GloMax Veritas Luminometer (Promega). 11 To calculate the proportion of viable cells, RLUs from the CTG assay were averaged between 12 replicates and normalized to the DMSO control. The proportion of inhibited cells was calculated as 13 one minus the proportion of viable cells. Drug synergy was determined by calculating the Bliss 14 excess score (Bliss 1956) 76 , i.e.

18
, Cleaved caspase 3/7 assay 21 Cells were seeded into 24-well plates at 100,000 cells/mL in triplicate and were treated with 22 decitabine (50 nM, 100 nM or 200 nM on days 0, 1 and 2) with and without RG3039 (2 uM on day 23 0). Cells were harvested on day 5 and the proportion of apoptotic cells was assessed using the 24 NucView 488 Caspase-3 Assay Kit (Biotium) according to the manufacturer's instructions and an 25 Attune NxT flow cytometer (Thermo Fisher Scientific) gating on the BL-1 channel.

27
Cell cycle assay 28 Cells were seeded into 24-well plates at 100,000 cells/mL in triplicate and were treated with 29 decitabine (50 nM, 100 nM or 200 nM on days 0, 1 and 2) with and without RG3039 (2 uM on day 30 0). Cells (500,000-1,000,000 per sample) were harvested on day 5 and the proportions of cells in 31 each phase of the cell cycle were assessed using the FxCycle Violet Kit (Thermo Fisher Scientific) 32 and an Attune NxT flow cytometer (Thermo Fisher Scientific) gating on the VL-1 channel. Briefly, 33 cells were washed once with PBS, fixed with 70% ethanol overnight at -20 ° C, pelleted, and then 34 washed with PBS 1-2 times. Cells were then resuspended in 1 mL permeabilization buffer (PBS 35 with 1% FBS and 0.1% Triton X-100) and 1 uL Fx cycle and stained for 30 minutes in the dark 36 before being analyzed via flow cytometry. 37 38 39 RNA-seq experimental procedures 1 2 3' RNA-seq 3 3' RNA-seq was performed to assess differential gene expression following decitabine and/or 4 RG3039 treatment. Cells were seeded into 6-well plates at 100,000 cells/mL in duplicate and were 5 treated with decitabine (100 nM on days 0, 1 and 2), RG3039 (2 uM on day 0), both drugs or DMSO. 6 On day 3, RNA was extracted using the RNeasy Mini Kit (Qiagen) according to the manufacturer's 7 instructions. RNA-seq libraries were prepared using the QuantSeq 3′ mRNA-Seq Library Prep Kit 8 FWD for Illumina (Lexogen) and assessed on a BioAnalyzer 2100 (Agilent) for library 9 quantification and quality control. RNA-seq libraries were sequenced on an Illumina HiSeq 4000 10 using single-end, 50-base pair sequencing.

12
Stranded RNA-seq 13 Stranded RNA-seq was performed for experiments in which strand directionality was required for 14 downstream analysis. Cells were seeded into 6-well plates at 100,000 cells/mL in duplicate or 15 triplicate and were treated with decitabine (100 nM on days 0, 1 and 2) or DMSO. The Salmon-tximport-DESeq2 pipeline 29 We used a workflow hereafter referred to as the "Salmon-tximport-DESeq2 pipeline" to perform 30 differential gene expression analysis. Salmon (version 1.2.1) was first used to quantify transcript 31 abundance 58 . A Salmon index was generated using the GENCODE (version 34) genome annotation, 32 and subsequently the `salmon quant` tool was used with the `--validateMappings` option to calculate 33 transcript abundances 77 . Then, the R package tximport was used to import Salmon results into R and 34 perform data preparation 57 . The `summarizeToGene` function was used to collapse transcript 35 abundances to the gene level. From here, the R package DESeq2 was used for differential gene 36 expression analysis 56 . We first extracted normalized counts for each RNA-seq experiment using 37 DESeq2 by running the `estimateSizeFactors` function and then the `counts` function with option 38 `normalized=TRUE`. For each individual experiment, the DESeq2 statistical model was modified 39 based on the experimental design. For experimental designs with multiple variables (e.g., multiple 40 drug conditions, time points, etc.), we used the likelihood ratio test (LRT) to perform differential 41 expression analysis. The LRT is conceptually similar to an analysis of variance (ANOVA) 42 calculation in a linear regression model 78 . In these cases, we specified the model design in the 1 `DESeq2` function as `~0 + variable1 + variable2 + variable1:variable2` and the argument 2 `test=LRT`. In simple experimental designs with one variable (e.g., DMSO vs. decitabine treatment), 3 DESeq2 was used with default options (i.e., a Wald test was used instead of a LRT). In these cases, 4 the model design was specified as `~cond`. For experiments with batch effects, the model design was 5 specified as `~cond + reps`. 6 7 Differential RNA stability analysis 8 9 The STAR-featureCounts-REMBRANDTS-limma pipeline 10 For analyses which required measurements of pre-mRNA and mature mRNA abundances from 11 RNA-seq samples (i.e., differential RNA stability analysis), we used a workflow hereafter referred to 12 as the "STAR-featureCounts-REMBRANDTS-limma pipeline". RNA-seq sequencing reads were 13 first aligned to the hg38 reference genome using STAR (version 2.7.3a) 59 . Then, featureCounts was 14 used to quantify intron and exon level counts. Finally, REMBRANDTS was used to calculate mRNA 15 stability as previously described (https://github.com/csglab/REMBRANDTS) 61 . Briefly, the package 16 estimates a gene-specific bias function that is subtracted from Δexon-Δintron calculations to provide 17 unbiased mRNA stability measurements. To assess differential RNA stability changes, we used 18 limma, which was designed for microarray experiments and serves a similar function to DESeq2, 19 though it supports negative values (relevant for RNA stability analysis) 60 . The model designs used 20 here are analogous to the designs for differential expression analysis described above. 21 22 RNA-seq iPAGE analysis 23 24 iPAGE run and downstream analysis pipeline 25 The iPAGE algorithm was used for gene set and pathway enrichment analysis on differential RNA 26 expression and stability results 62 . MSigDB (version 7.4.) was downloaded and modified to be 27 compatible with iPAGE workflow 33,71,72 . iPAGE was used in continuous mode, which accepts gene-28 level numeric values (e.g., logFCs) as input. Briefly, iPAGE quantizes differential measurements 29 into equally populated bins and then, for every given geneset, calculates the mutual information (MI) 30 between each cluster bin and a binary vector of pathway memberships for genes in a given gene set.

31
The significance of each MI value is then assessed through a randomization-based statistical test and 32 hypergeometric distribution to determine whether there is over or under representation of a gene set 33 in each cluster bin. Pre-processing HERV annotations for alignment tasks 1 2 Annotations in BED12 format were downloaded from the Human Endogenous RetroViruses 3 Database 45,46 . To prepare these annotations for alignment tasks, i.e. building Salmon and STAR 4 indices, CGAT Apps was used to convert BED12 files to GTF format (`cgat bed2gff --as-gtf`) and 5 the `getfasta` module from bedtools (with options `-name+ -split`) was used to convert BED12 files 6 to FASTA format 79-80 . Reproducible scripts for preparing ERV annotations for alignment tasks are 7 available at https://github.com/abearab/HERVs. 8 9 RNA-seq workflows for specific experiments 10 11 Decitabine and RG3039 drug combination experiments 12 We performed 3' RNA-seq on HL-60 cells treated with DMSO, decitabine alone, RG3039 alone or 13 both drugs for 72 hours in duplicate (see above for experimental procedures). Raw sequencing data 14 were processed using our Salmon-tximport-DESeq2 pipeline (see above). DESeq2 was used to 15 conduct differential gene expression analysis using a likelihood ratio test and the model design `~0 + 16 decitabine + rg3039 + decitabine:rg3039`. Pathway enrichment was assessed using iPAGE (see 17 above). For PCA analysis, the `varianceStabilizingTransformation` function from the DESeq2 18 package was used to prepare counts. The `plotPCA` function was used to calculate PC variances as 19 percentages. Finally, `ggplot2` was used to visualize a two-dimensional representation of the PCA 20 analysis. For differential ERV expression analysis, processed ERV annotations (see above) in 21 FASTA format used to build an index for Salmon workflow and then samples processed through the 22 Salmon-tximport-DESeq2 pipeline (see above). Upregulated ERVs were defined as p-value < 0.05 23 and log2FC > 2, and downregulated ERVs were defined as p-value < 0.05 and log2FC < -2. The 24 intersections of ERV data were visualized using UpSet plots in Python 81 . 25 26

HL-60 time-series experiments 27
We performed stranded RNA-seq on HL-60 cells treated with decitabine for 6 hours, 72 hours and 28 120 hours in duplicate (see above for experimental procedures). Differential expression analysis was 29 performed using our Salmon-tximport-DESeq2 pipeline (see above), using a likelihood ratio test and 30 a two-variable model design incorporating treatment condition (decitabine or DMSO) and time (6, 31 72 or 120 hours). Differential RNA stability analysis was performed using our STAR-featureCounts-32 REMBRANDTS-limma pipeline (see above). Pathway enrichment for differential expression and 33 RNA stability data was assessed using iPAGE (see above).

35 AML cell line panel experiments 36
We performed stranded RNA-seq on AML cell lines treated with decitabine or DMSO for 72 hours 37 in three replicates (see above for experimental procedures). Differential expression analysis was 38 performed using our Salmon-tximport-DESeq2 pipeline (see above), using a Wald test. Differential 39 RNA stability analysis was conducted using our STAR-featureCounts-REMBRANDTS-limma 40 pipeline (see above). Pearson correlation tests from the Hmisc and corrplot R packages were used to 41 assess correlation between differentially expressed genes in HL-60 and other AML cell lines. UpSet 42 plots in Python were used to identify and visualize genes across multiple cell lines that conferred 1 drug sensitivity in the CRISPRi screen (ρ score < -0.1 and p < 0.05), were RNA hypermethylated 2 (log2FC > 1 and p < 0.05) upon decitabine treatment, and either had decreased expression or RNA 3 stability (log2FC < -0.1 and p < 0.05) upon decitabine treatment 81 . 4 5 MeRIP-seq 6 7 Experimental procedure 8 We performed MeRIP-seq as previously described on HL-60 cells treated with DMSO or decitabine 9 for 72 hours in biological duplicates 37 . First, 2 µg of the fragmented total RNA per sample was used 10 for RNA immunoprecipitation (IP) with 5 µg of the anti-m 6 A antibody (ABE572, Millipore). RNA-11 seq libraries from input and IP samples were prepared using the SMARTer Pico Input Mammalian 12 v2 RNA-seq kit (Takara) and sequenced as SE50 runs on an Illumina HiSeq4000.

14
Alignment task for mRNAs of protein coding genes and ERVs 15 MeRIP-seq reads were aligned to the hg38 reference genome using STAR ( options `--method=merge --merge-by-name`) and the `getfasta` module from bedtools (with options 27 `-name -s -split`) were used for sequence extraction 79,80 . Finally, the FIRE algorithm was used in 28 non-discovery mode for enrichment analysis of known m 6 A motifs (i.e., RGAC or [AG]GAC) within 29 peak sequences, compared to randomly generated sequences 39 .

31
Peak calling and differential RNA methylation analysis 32 RADAR was used to perform peak calling and differential methylation analysis 41 . Differentially 33 methylated peaks were defined as FDR < 0.1 and logFC > 0.5. The logFC values for protein coding 34 genes and each of ERVs used to test global hypermethylation using Wilcoxon test and t-test 35 functions with `mu=0`, `alternative="greater"` options. Results are shown as annotated volcano plots 36 using ggplot2 in R. For peak visualization across individual mRNA transcripts, the `plotGeneCov` 37 function from the RADAR R package was used to generate coverage plots. Then, the Gviz R 38 Bioconductor package was used to draw detailed information for each mRNA transcript 84 .

12
Monosomes were isolated using MicroSpin S-400 HR (Cytiva) columns, pre-equilibrated with 3 mL 13 polysome buffer per column. 100 µl digested lysate was loaded per column (two columns were used 14 per 200 µl sample) and centrifuged 2 min at 600 g. The RNA from the flow through was isolated 15 using the RNA Clean and Concentrator-25 kit (Zymo). In parallel, total RNA from undigested 16 lysates were isolated using the same kit.

18
Ribosome protected footprints (RPFs) were gel-purified from 15% TBE-Urea gels as 17-35 nt 19 fragments. RPFs were then end-repaired using T4 PNK (New England Biosciences) and pre-20 adenylated barcoded linkers were ligated to the RPFs using T4 Rnl2(tr) K227Q (New England 21 Biosciences). Unligated linkers were removed from the reaction by yeast 5'-deadenylase (New 22 England Biosciences) and RecJ nuclease (New England Biosciences) treatment. RPFs ligated to 23 barcoded linkers were pooled, and rRNA-depletion was performed using riboPOOLs (siTOOLs) per 24 the manufacturer's recommendations. Linker-ligated RPFs were reverse transcribed with ProtoScript 25 II RT (New England Biosciences) and gel-purified from 15% TBE-Urea gels. cDNA was then 26 circularized with CircLigase II (Epicentre) and used for library PCR. First, a small-scale library PCR 27 was run supplemented with 1X SYBR Green and 1X ROX (Thermo) in a qPCR instrument. Then, a 28 larger scale library PCR was run in a conventional PCR instrument, performing a number of cycles 29 that resulted in ½ maximum signal intensity during qPCR. Library PCR was gel-purified from 8% 30 TBE gels and sequenced on a SE50 run on an Illumina HiSeq4000.

32
Data preprocessing 33 The adapters in the sequencing reads were removed using cutadapt 86 (v3. Ribolog was used to compare translational efficiency across conditions 7 (https://github.com/goodarzilab/Ribolog) 92 . Briefly, Ribolog applies a logistic regression to model 8 individual Ribo-seq and RNA-seq reads in order to provide estimates of logTER (i.e., logFC in TE) 9 and its associated p-value across the coding transcriptome. 10 11 Multiomics data integration 12 13 To identify candidate genes among our multiomics datasets for downstream validation of our 14 decitabine-m 6 A model, we examined the intersection of three sets of genes: (1) sensitizing hits in the 15 CRISPRi screen, defined as ρ score < -0.1 and p < 0.05; (2) genes with downregulated expression 16 upon decitabine treatment, defined as log2FC < -0.1 and p < 0.05; (3) genes with RNA 17 hypermethylation upon decitabine treatment, defined as logFC > 1 and p < 0.05. Intersections 18 between sets were visualized through a Venn diagram in Python.

20
Quantitative RT-PCR 21 22 RNA isolation 23 Total RNA for quantitative RT-PCR assays was isolated using the Quick-RNA Microprep kit 24 (Zymo) with on-column DNase treatment per the manufacturer's protocol.

26
Quantitative RT-PCR 27 Transcript levels were measured using quantitative RT-PCR by reverse transcribing total RNA to 28 cDNA (Maxima H Minus RT, Thermo Fisher Scientific), then using fast SYBR green master mix 29 (Applied Biosystems) or Perfecta SYBR green supermix (QuantaBio) per the manufacturer's 30 instructions. HPRT1 was used as an endogenous control. 31 32

11
(G) A cleaved caspase 3/7 assay shows the fraction of apoptotic HL-60 cells at day 5 following treatment with

12
DMSO or decitabine ± RG3039. Data are shown as means ± SD for three replicates.       to capture abundances of pre-mRNA (left), mature mRNA (middle) and predicted m 6 A hypermethylated loci 14 for each gene (right). Data are shown as three replicates, and one-tailed Mann-Whitney U-tests were used to 15 assess statistical significance.

13
Data are shown as means ± SD for three replicates. Data were derived from the same experiment as Figure   14 1G.

16
A B

Retroposons SINEs LINEs
Differential RNA methylation