{'fraction_non_alphanumeric': 0.04344005904474045, 'fraction_numerical': 0.03124450848768144, 'mean_word_length': 4.8922137088424105, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 2, 'https://': 10, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
N 6 -methyladenine DNA Modification in Glioblastoma

Genetic drivers of cancer can be dysregulated through epigenetic modifications of DNA. Although the critical role of DNA 5-methylcytosine (5mC) in the regulation of transcription is recognized, the functions of other non-canonical DNA modifications remain obscure. Here, we report the identification of novel N6-methyladenine (N6-mA) DNA modifications in human tissues and implicate this epigenetic mark in human disease, specifically the highly malignant brain cancer glioblastoma. Glioblastoma markedly upregulated N6-mA levels, which co-localized with heterochromatic histone modifications, predominantly H3K9me3. N6-mA levels were dynamically regulated by the DNA demethylase ALKBH1, depletion of which led to transcriptional silencing of oncogenic pathways through decreasing chromatin accessibility. Targeting the N6-mA regulator ALKBH1 in patient-derived human glioblastoma models inhibited tumor cell proliferation and extended the survival of tumor-bearing mice, supporting this novel DNA modification as a potential therapeutic target for glioblastoma. Collectively, our results uncover a novel epigenetic node in cancer through the DNA modification N6-mA.


INTRODUCTION
Aberrant epigenetic landscapes promote tumor initiation and progression (Elsä sser et al., 2011;Lové n et al., 2013;Northcott et al., 2014). To date, the focus of cancer research has been on the global and local aberrations of 5-methylcytosine (5mC). Genome-wide 5mC hypomethylation frequently occurs in cancer genomes, leading to widespread genomic instability and derepression of repetitive elements (Cadieux et al., 2006;Feinberg and Vogelstein, 1983;Widschwendter et al., 2004). However, local hypermethylation of 5mC events is often found in close proximity to important tumor suppressor genes in a number of cancers (Esteller, 2002). Low-grade gliomas and secondary glioblastomas commonly harbor mutations in isocitrate dehydrogenase 1 or 2 (IDH1/2), leading to DNA hypermethylation and manifesting as the glioma-CpG island methylator phenotype (G-CIMP) (Noushmehr et al., 2010;Turcan et al., 2012). Thus, DNA methylation is a critical epigenetic mark that regulates many important developmental processes and plays a fundamental role in cancer.
Although numerous epigenetic modifications have been described in addition to DNA methylation, it is still unknown whether tumor cells can adopt novel epigenetic mechanisms that are rarely utilized in normal adult tissues. In addition to the canonical 5mC, other DNA methylation events, including N 6 -methyladenine (N 6 -mA), have been identified in bacteria (Vanyushin et al., 1968) and a limited number of eukaryotes, such as mosquitoes, plants, C. elegans (A) Levels of the N 6 -methyladenine (N 6 -mA) DNA modification were assessed via DNA dot blot in normal human astrocytes, patient-derived GSC models (387, D456, GSC23, and 1919), and primary human glioblastoma specimens (3028 and CW2386) using an N 6 -mA-specific antibody. Methylene blue detected DNA loading. (B) Mass spectrometry analysis of N 6 -mA in two normal human astrocyte cell lines and two patient-derived GSC models (387 and D456). Data are presented as mean ± SD. Two replicates were used for each sample. Significance was determined by one-way ANOVA with Tukey multiple comparisons test. p < 0.0001 for each human astrocyte versus GSC comparison. (C) N 6 -mA DNA immunofluorescence in normal human astrocytes and human patient-derived GSC models (387, D456, and GSC23). DAPI indicates cell nuclei. Scale bars, 50 mm. et al., 2017;Zhang et al., 2015). We recently demonstrated that N 6 -mA DNA methylation occurs in mammals, where it acts mainly as a repressive mark to silence transcription of long interspersed nuclear element (LINE) transposons in mouse embryonic stem cells (ESCs) (Wu et al., 2016), and this modification has also been reported in human tissue (Xiao et al., 2018;Zhu et al., 2018). However, the roles of N 6 -mA in human cancer and in genome biology are largely unknown. To interrogate N 6 -mA function, we utilized glioblastoma as a model human cancer. Glioblastoma is the most prevalent and aggressive primary intrinsic brain tumor and is characterized by widespread epigenetic dysregulation, including DNA hypermethylation of 5mC and alterations in chromatin remodeling enzymes, such as the polycomb complex proteins EZH2 and BMI1 (Abdouh et al., 2009;Jin et al., 2017;Suvà et al., 2009). Thus, we studied the N 6 -mA DNA modification in human patient-derived glioblastoma tumor models and functionally characterized its genomic localization and role in the regulation of glioblastoma.

Identification of N 6 -mA DNA Modifications in Human Glioblastoma
To investigate the roles of the N 6 -mA DNA modification in human glioblastoma, we extracted genomic DNA from functionally validated human patient-derived glioblastoma stem cell (GSC) models (387, D456, GSC23, and 1919) and primary human tumor specimens (3028 and CW2386) and then performed dot blot analysis using an N 6 -mA-specific antibody that has been validated previously (Greer et al., 2015;Wu et al., 2016;Zhang et al., 2015). GSCs and primary human tumor specimens displayed strikingly elevated N 6 -mA levels relative to normal human astrocytes (Figure 1A). N 6 -mA levels were independently quantified by mass spectrometry (MS) using an established and highly sensitive (limit of quantification [LOQ], 1.6 fmol) MS (ultra-high performance liquid chromatography tandem-mass spectrometry [UHPLC-MS/MS]) approach with stable isotope-labeled [ 15 N 5 ] N 6 -mA as an internal standard for sample enrichment and quantification (Figures S1A and S1B), which demonstrated that N 6 -mA levels in GSCs were elevated more than 100-fold compared with normal human astrocytes (Figure 1B). N 6 -mA DNA levels were also measured by immunofluorescence staining. To eliminate signal from RNA, samples were treated with RNaseA prior to incubation with primary antibodies. The N 6 -mA DNA levels detected by immunofluorescence validated the dot blot and MS results, with much higher N 6 -mA levels in patient-derived GSCs than in normal human astrocytes ( Figures 1C and 1D). Effects of cell culture conditions were ruled out using DNA N 6 -mA immunohistochemistry on a tissue microarray (TMA) with normal brain and glioblas-toma tissues, confirming increased DNA N 6 -mA levels in glioblastoma (Figures 1E and 1F). Potential inter-individual heterogeneity in N 6 -mA levels was addressed using immunofluorescence staining of matched glioblastoma and surrounding normal brain tissue from the same patient following surgical resection. Consistently, glioblastoma tumor tissue contained increased levels of N 6 -mA (Figure S1C). N 6 -mA levels were elevated in other CNS cancers, including diffuse intrinsic pontine glioma (DIPG), meningioma, and medulloblastoma, compared with normal human astrocytes ( Figure S1D). Of note, a recent study reported high levels of N 6 -mA in normal human adult tissues (at several hundred parts per million), with downregulation in gastric and stomach cancers (Xiao et al., 2018). However, our current study, together with other previous studies , did not find high levels of N 6 -mA in normal adult tissues or mammalian cells. Collectively, our results support elevation of a novel DNA modification, N 6 -mA, in glioblastoma, with likely involvement in other cancers.
Genomic Profiling and Pattern of N 6 -mA in Glioblastoma N 6 -mA genomic localization was analyzed using N 6 -mA DIP sequencing (DIP-seq) (DNA immunoprecipitation with anti-N 6 -mA antibodies followed by next-generation sequencing) in two human patient-derived GSC models (387 and GSC23) and one primary human tumor specimen (CW2386) to identify N 6 -mA-enriched genomic regions. The number of N 6 -mA peaks measured by DIP-seq in glioblastoma ranged from 7,282 to 17,263 per model (Figure 2A). N 6 -mA peaks were most common in intergenic regions, consistent with our previous findings in murine ESCs (Wu et al., 2016;Figure 2B). N 6 -mA peaks were present on each chromosome, with mild enrichment on chromosomes 7, 19, and 21 (Figure 2C). Gene ontology (GO) enrichment analysis of the closest gene to each of the shared N 6 -mA peaks revealed an enrichment in neurogenesis and neuronal development pathways (Figures 2D and 2E). GO analysis restricted to peaks within gene bodies revealed similar results (Figures S2A and S2B). Although adenine nucleotides are distributed throughout the genome, selective distribution of N 6 -mA could reflect sequence-specific recognition of molecular regulators. To investigate this, we identified DNA motifs enriched within N 6 -mA regions in the common shared peaks. One of the topranking motifs (GGAAT) closely resembles human major satellite repeats that are enriched at constitutive heterochromatic regions (Figure S2C). Consistent with findings that DNA N 6 -mA functions as a repressive mark in murine ESCs and brains (Wu et al., 2016;Yao et al., 2017), these data suggest that N 6 -mA-mediated repression of neurodevelopmental pathways may serve a role in heterochromatin formation that contributes to glioblastoma tumorigenesis. (D) Quantification of the percentage of N 6 -mA-positive cells by immunofluorescence staining in (C). N 6 -mA was quantified by counting 100 cells from each sample. n = 3 slides/cell type. Data are presented as mean ± SD. Significance was assessed by ANOVA. ***p < 0.001. (E) Immunohistochemistry (IHC) staining of N 6 -mA in non-neoplastic brain tissues (6 total) and human primary glioblastoma specimens (67 total) from a tissue microarray. Scale bars, 50 mm. (F) Quantification of N 6 -mA levels from the immunohistochemistry staining in (E). N 6 -mA levels were scored from low levels (score = 0) to highest levels (score = 3). Data are presented as mean ± SD. Student's t test, p = 0.0006. See also Figure S1.   Figure 2. Genomic Localization of N 6 -mA Enrichment in Glioma Stem Cells (A) N 6 -mA enrichment was identified using the spatial clustering for identification of ChIP-enriched regions (SICER) algorithm for a human primary glioma sample (CW2386) and two in vitro human patient-derived glioma models (T387 and GSC23). The number of N 6 -mA peaks are shown, with the number of shared peaks in the center. (B) Genome ontology analysis showing the fraction of common N 6 -mA peaks present in distal intergenic, intronic, or gene regions compared with the genome background (chi-square test, p < 0.0001). (C) Genome ontology analysis showing the fraction of common N 6 -mA peaks on each chromosome compared with the genome background. Chromosome 7, p = 1.2 3 10 À208 ; chromosome 21, p = 2.2 3 10 À322 ; chromosome 3, p = 5.1 3 10 À68 ; chromosome 5, p = 1.2 3 10 À69 . (D) Gene ontology (GO) enrichment analysis of the closest gene to each of the common N 6 -mA peaks. Values are expressed as -log10 (false discovery rate [FDR]corrected p value). (E) Enrichment map demonstrating key pathways identified in the GO enrichment analysis in (D). Circles represent individual gene sets. The size of the circle depicts the number of genes in the gene set, and the edge color depicts the FDR-corrected p value of the enrichment, with dark blue representing the most significant gene sets. See also Figure S2.  3B). Aggregation analysis showed that N 6 -mA peaks common to all three glioma models ( Figure 2A) were directly superimposed with those of the heterochromatic histone marks ( Figure 3C). All N 6 -mA peaks in the 387 patient-derived GSC model, some of which were unique in this model, were strongly correlated with these heterochromatic marks ( Figure 3D). Genomic overlap between the N 6 -mA DNA modification and the canonical 5mC DNA methylation mark was interrogated through whole-genome bisulfite sequencing analysis of the 387 patient-derived GSCs and mapping of the intersection with N 6 -mA. As expected, 5mC marks were depleted in promoter regions but slightly increased in N 6 -mA peaks compared with randomly selected genomic regions ( Figures S2D-S2G). Thus, our data suggest that N 6 -mA strongly co-localizes with heterochromatic histone modifications and with the largely repressive canonical DNA methylation mark 5mC. Because the N 6 -mA DNA modification acts as a repressive mark that is highly associated with the H3K9me3 heterochromatin histone modification, the high levels of N 6 -mA may function to support GSC survival by repressing key tumor suppressor genes. To investigate this possibility, the Tumor Suppressor Gene (TSG) Database (Zhao et al., 2016) was leveraged to identify potential tumor suppressor genes that contain both N 6 -mA and H3K9me3 enrichment. Several tumor suppressor genes, including CDKN3, RASSF2, and AKAP6, appear to be repressed by N 6 -mA and H3K9me3 ( Figure S3). ALKBH1 Is a DNA N 6 -mA Demethylase that Actively Modulates Gene Expression in Human Glioblastoma DNA methylation is dynamically regulated by a number of enzymes. We next investigated the enzymes controlling DNA N 6 -mA modification in glioblastoma. Our previous study in murine ESCs identified alkylated DNA repair protein AlkB homolog (ALKBH1) as a demethylase for DNA N 6 -mA (Wu et al., 2016), so we focused on its human homolog as a potential regulatory node. In vitro DNA demethylation assays demonstrated strong N 6 -mA demethylase catalytic activity of the human ALKBH1 protein toward N 6 -mA in oligonucleotides using dot blotting (Figure 4A) or UHPLC-MS ( Figure 4B). To investigate its functions in glioblastoma, ALKBH1 expression was targeted by two independent, non-overlapping small hairpin RNA (shRNA) lentiviral constructs (designated shALKBH1.1008 and shALKBH1.1551) compared with a non-targeting control sequence shRNA insert (shCONT), which increased DNA N 6 -mA levels in two patientderived GSCs (387 and D456; Figure 4C). Because a recent study indicated that ALKBH1 might have activities against tRNA (Liu et al., 2016), we further investigated and found that ALKBH1 does not function as a tRNA-1mA demethylase in patient-derived GSCs (data not shown).
We next aimed to identify N 6 -mA methyltransferases in GSCs. N-6 adenine-specific DNA methyltransferase 1 (N6AMT1) was recently suggested as a DNA N 6 -mA methyltransferase in human cells (Xiao et al., 2018). In light of these findings, we interrogated the function of N6AMT1 in patient-derived GSCs. In contrast to this previous report, genomic knockout of N6AMT1 in GSCs did not alter the levels of DNA N 6 -mA, as assessed by dot blot (Figures S4A and S4B). Further, purified N6AMT1 recombinant protein incubated with DNA oligonucleotides in vitro did not alter DNA N 6 -mA levels, suggesting that N6AMT1 does not serve as an N 6 -mA methyltransferase, at least in glioblastoma ( Figures S4C and S4D). Our findings are consistent with a previous study that also did not support N6AMT1 as an N 6 -mA methyltransferase (Liu et al., 2010). Thus, we focused our attention on ALKBH1, the only definitive modulator of N 6 -mA in GSCs, to understand the role of N 6 -mA regulation in glioblastoma.
Our data suggest that ALKBH1 may control specific genetic pathways in glioblastomas. To investigate this, we performed RNA sequencing (RNA-seq) to interrogate the effects of ALKBH1 knockdown on the whole transcriptome. 321 genes were downregulated by 2-fold or greater, whereas only 37 genes were upregulated by more than 2-fold upon ALKBH1 depletion (Figure 4D), suggesting that an increase in N 6 -mA DNA levels primarily represses gene expression in human glioblastoma. To further elucidate the genomic sites of ALKBH1 demethylase activity, N 6 -mA DIP-seq was performed on patient-derived GSCs (387) after transduction with either one of two independent ALKBH1 shRNAs or shCONT. 68% of downregulated genes  (A) The fraction of N 6 -mA peaks in a patient-derived GSC model (387) that overlap with matched H3K9me3, H3K27me3, and H3K4me3 histone modification domains in the same model. (B) Number of overlaps between N 6 -mA peaks and histone modification domains in a patient-derived GSC model (387) relative to the genome background. (C) Heatmap and profile plot demonstrating the intersection of N 6 -mA common peaks with H3K9me3 and H3K27me3 signals over a scaled window 5 kb upstream and downstream of the common N 6 -mA peak. (D) Heatmap of N 6 -mA peaks in a GSC model (387) divided into those that are cobound with H3K9me3 and H3K27me3, those bound by H3K9me3 or H3K27me3 alone, and those bound by H3K4me3 alone. The signal is shown over a scaled window 5 kb upstream and downstream of the N 6 -mA peak, and the height of the heatmap is directly proportional to the number of regions present in each segment. See also Figure S3.  after ALKBH1 knockdown were associated with genomic regions with increased N 6 -mA levels ( Figure 4E), strongly suggesting that the transcriptional downregulation is a direct consequence of elevated N 6 -mA levels because of ALKBH1 depletion. To understand the role of N 6 -mA in regulating chromatin accessibility, ATAC-seq (assay for transposase-accessible chromatin using sequencing) was performed following knockdown of ALKBH1. Consistent with our finding of co-localization of N 6 -mA with heterochromatic regions, ALKBH1 knockdown led to a greater number of genomic loci with decreased chromatin accessibility (1,389 sites) than sites with increased accessibility (243 sites) ( Figures 4F-4I). ATAC-seq peaks, representing sites of open chromatin, negatively correlated with N 6 -mA peaks, indicating that N 6 -mA is depleted at sites of high chromatin accessibility (Figures 4J). Overlayed ATAC-seq and RNA-seq data showed that more genes bearing a significantly decreased accessibility site within 2 kb of their transcriptional start sites tended to be transcriptionally repressed (37 genes) rather than upregulated (19 genes) ( Figure 4K). To further investigate the interaction between ALKBH1 and N 6 -mA modifications, biotin-labeled oligonucleotides with or without N 6 -mA modifications were used to perform ALKBH1 pull-downs using whole-cell lysates from patient-derived GSCs. N 6 -mA modified oligonucleotides preferentially bound ALKBH1 relative to unlabeled oligonucleotides, suggesting that ALKBH1 binds regions bearing N 6 -mA modifications (Figure S5A). To assess the preference of ALKBH1 to bind N 6 -mA genome-wide, ALKBH1 ChIP-seq was used to identify ALKBH1 binding sites. ALKBH1 ChIP-seq peaks were preferentially enriched within N 6 -mA marked regions compared with the genome background ( Figures S5B and S5C). Thus, ALKBH1 functions as a transcriptional activator by binding to regions of N 6 -mA enrichment and removing repressive N 6 -mA marks at selected genomic loci.

ALKBH1 Depletion Facilitates Heterochromatin
Formation in Human Glioblastoma through N 6 -mA DNA Modification We next analyzed co-localization of N 6 -mA peaks and genomic sites directly regulated by ALKBH1 (i.e., N 6 -mA peaks gained after ALKBH1 knockdown). Consistent with the strong association of N 6 -mA and heterochromatin marks at a genome-wide level, N 6 -mA sites directly increased upon ALKBH1 knockdown were even more strongly correlated with heterochromatin histone modifications, particularly H3K9me3 ( Figures 5A and 5B). Genes marked by N 6 -mA peaks overlapped with those marked by H3K9me3 ( Figures 5C, 5D, and S5D). For example, two genes downregulated upon ALKBH1 depletion, RFTN2 and LOX, were marked by increased N 6 -mA levels that co-localized with H3K9me3-demarcated heterochromatin regions (Figures S5G and S5H).
To assess the effects of ALKBH1-mediated alterations of N 6 -mA on the global H3K9me3-marked heterochromatin landscape, we performed H3K9me3 ChIP-seq following ALKBH1 knockdown in patient-derived GSCs, revealing an increased number of H3K9me3 peaks following ALKBH1 knockdown (9,999 peaks with increased H3K9me3 signal versus 2,995 peaks with decreased H3K9me3 signal) ( Figures 5E and 5F), which was confirmed in two additional patient-derived GSCs ( Figures S5E and S5F). The genes SDK1 and BMT2 gained N 6 -mA peaks co-localizing with gained H3K9me3 signal following ALKBH1 knockdown ( Figure 5G). Collectively, these results suggest that N 6 -mA plays an important role in facilitating heterochromatin formation in human glioblastoma, particularly acting through the H3K9me3 histone modification.

ALKBH1 Regulates Hypoxia Response Genes
To interrogate the gene pathways regulated by ALKBH1-sensitive N 6 -mA sites, we performed GO analysis on genes  (B) The in vitro demethylation reaction was quantified by liquid chromatography-tandem mass spectrometry (LC-MS/MS) following addition of ALKBH1 protein to N 6 -mA-labeled DNA oligonucleotides. Data are presented as mean ± SD (Student's t test, ***p < 0.001, n = 3). (C) ALKBH1 expression was decreased using two independent shRNAs in human patient-derived GSC models, and N 6 -mA levels were assessed using a DNA dot blot. Methylene blue detected DNA loading. The ALKBH1 protein level was assessed using western blot. (D) RNA-seq analysis in a human patient-derived GSC model following ALKBH1 knockdown. Blue dots indicate the most highly upregulated genes following ALKBH1 knockdown (37), whereas red dots indicate the most highly downregulated genes following ALKBH1 knockdown (321). (E) Venn diagram indicating overlap between (1) genomic regions with gained N 6 -mA after ALKBH1 knockdown and (2) downregulated genes after ALKBH1 knockdown. (F) Genome-wide differentially accessible sites were identified by ATAC-seq. 1,632 sites of differential accessibility (FDR-corrected p < 0.05) were identified and are visualized in an MA plot. 1,389 sites were identified with decreased accessibility after ALKBH1 knockdown, whereas 243 sites displayed increased accessibility. Log2 fold change > 0.5; p < 0.05. (G) Chromatin accessibility heatmap for the differentially accessible sites (absolute value of log2 fold change > 1). Two replicates were performed for each treatment group. The signal is shown over a scaled window 5 kb upstream and downstream of the differentially accessible region. (H) Supervised heatmap showing the correlation between the ATAC-seq counts in non-targeting and ALKBH1 knockdown samples. Two replicates were performed for each sample. (I) Profile plot showing the ATAC-seq signal over sites of decreased chromatin accessibility in ALKBH1 knockdown samples. The scaled signal is shown 5 kb upstream and downstream of each ATAC peak. (J) Heatmap and profile plot showing the N 6 -mA DIP-seq and ATAC-seq signals over the top 10,000 ranked ATAC-seq peaks in the 387 GSC model. (K) Graph showing the mRNA fold change following ALKBH1 knockdown for genes with sites of decreased accessibility. Among genes with decreased accessibility sites within 2 kb of their transcriptional start site, 37 genes were downregulated, and 19 genes were upregulated. See also Figure S4.  downregulated upon ALKBH1 depletion, revealing enrichment of key oncogenic pathways previously implicated in glioblastoma pathogenesis, including hypoxia responses, nervous system development, cellular and tissue development, and plasma membrane structure pathways (Figures S6A and S6B). Of note, hypoxia response genes were modulated by ALKBH1 even when cells were grown under room air conditions. Because hypoxia-related gene signatures were the top enrichment following ALKBH1 knockdown, the comparative transcriptional effects of ALKBH1 knockdown were assessed under hypoxic conditions (1% oxygen) by RNA-seq. Similar to the normoxic conditions, genes in the hypoxia pathway were downregulated following ALKBH1 knockdown when cells were grown in hypoxia (Figures S6C and S6D). Upregulated genes upon ALKBH1 knockdown in hypoxia included DNA damage and p53 pathway genes (Figures S6E and S6F). To support the clinical significance of these findings, an ALKBH1-regulated gene signature highly correlated with a hypoxia signature in The Cancer Genome Atlas (TCGA) glioblastoma dataset (Figure S6G). To explore the correlation between hypoxia and ALKBH1-dependent N 6 -mA demethylation sites, RNA-seq of GSCs under hypoxia identified that over 50% of genes upregulated in hypoxia were suppressed by ALKBH1 knockdown (Figure S6H). N 6 -mA peaks that decreased following hypoxia tended to be located closer to upregulated genes in hypoxia (e.g., MIAT) than randomly selected peaks (genome background) ( Figures S6I-S6K). Thus, ALKBH1 regulates key processes previously implicated in GSC biology through alterations in the N 6 -mA landscape. H3K9me3 has been previously linked to repression of hypoxia-induced genes during development and tumorigenesis (Luo et al., 2012), including in glioblastoma (Intlekofer et al., 2015). Collectively, our results demonstrate an unexpected role of N 6 -mA and ALKBH1 in modulating hypoxia-induced genes in human glioblastoma. N 6 -mA Is a Potential Therapeutic Target in Glioblastoma To address the functional significance of N 6 -mA in tumor biology, we investigated ALKBH1 as a potential therapeutic target. ALKBH1 protein levels were modestly but consistently elevated in GSCs compared with their matched differentiated non-stem tumor cells by immunoblotting ( Figure S7A). Silencing ALKBH1 decreased GSC proliferation ( Figure 6A). To control for the potential off-target effects of shRNAs in additional experi-ments, we depleted ALKBH1 using single guide RNAs (sgRNAs) by CRISPR-Cas9 technology with two independent non-overlapping sgRNAs (sgALKBH1#4 and sgALKBH1#5). Following ALKBH1 depletion with sgRNAs, GSCs displayed reduced cell growth, consistent with the shRNA results ( Figure 6B). Sphere formation is surrogate marker of self-renewal, albeit with caveats. In vitro limiting dilution sphere formation assays revealed that ALKBH1 knockdown resulted in a more than 10-fold decrease in the frequency of sphere formation ( Figure 6C). GSCs are functionally defined by their ability to initiate tumors in vivo, prompting us to evaluate the potential anti-tumor effects of ALKBH1 targeting in vivo. Two different patient-derived GSCs (387 and D456) transduced with either one of two non-overlapping ALKBH1-targeting shRNAs (shALKBH1) or shCONT were implanted into the brains of immunocompromised mice. Consistent with our in vitro results, animals bearing tumors derived from GSCs expressing shRNAs or sgRNAs targeting ALKBH1 displayed increased survival relative to those bearing GSCs expressing shCONT or control sgRNA (sgCONT), respectively ( Figures 6D and 6E). In light of the role of ALKBH1 in the regulation of hypoxia pathways, we determined the effects of ALKBH1 knockdown under hypoxic conditions. Depleting ALKBH1 mRNA levels by shRNAs in hypoxia revealed reduced cell viability in vitro, similar to under normoxia ( Figure S7B). ALKBH1 knockdown in hypoxia impaired sphere formation capacity in a limiting dilution assay, mirroring our findings under normoxia ( Figure S7C). Overexpression of ALKBH1 diminished N 6 -mA levels in GSCs ( Figure S7D). In contrast to the knockdown studies, ALKBH1 overexpression did not alter cell viability, self-renewal, or in vivo tumor formation ( Figures S7E-S7G). Although N 6 -mA was reduced initially, its levels returned to baseline after 3 days, despite continuous overexpression of ALKBH1 (data not shown). This implies that compensatory mechanisms, which likely include currently unknown DNA demethylases and methyltransferases, exist to maintain the levels of N 6 -mA within a relatively narrow range that is compatible with cellular survival. In summary, ALKBH1 is essential for maintaining cell viability and stemness properties under both normoxic and hypoxic conditions.

17%
To determine the clinical relevance of these findings, we performed in silico studies on the TCGA glioblastoma dataset. Although correlation with patient outcome does not universally indicate the importance of any individual gene targets (Kaelin, 2017), ALKBH1 was highly expressed in glioblastomas relative Figure 5. ALKBH1 Regulates Downstream Gene Expression through N 6 -mA-Dependent Heterochromatin Formation (A) Pie chart depicting the fraction of N 6 -mA peaks gained after ALKBH1 knockdown in a patient-derived GSC model (387) that overlap with matched H3K9me3, H3K27me3, and H3K4me3 histone modification domains in the same model. (B) Overlaps between N 6 -mA peaks gained following ALKBH1 knockdown and histone modification domains in a patient-derived GSC model (387) relative to the genome background; **p < 0.01. (C) Pie chart showing the percentage of genes targeted by the N 6 -mA DNA modification that overlap with the H3K9me3 histone modification. (D) Heatmap showing gained N 6 -mA peaks after ALKBH1 knockdown with shRNA (shALKBH1.1008) and co-localization with H3K9me3. The signal is shown over a scaled window 10 kb upstream and downstream of the gained N 6 -mA peak. (E) Volcano plot showing the number of H3K9me3 peaks gained (red, 9,999) and lost (blue, 2,995) following ALKBH1 knockdown in the 387 GSC model. (F) Heatmap indicating sites of H3K9me3 enrichment following ALKBH1 knockdown with shALKBH1.1551. The signal is shown over a scaled window 5 kb upstream and downstream of the gained N 6 -mA peak.  Figure 6. ALKBH1 Is Essential for Glioblastoma Stem Cell Growth, Self-Renewal, and Tumor Formation Capacity (A) Top: cell proliferation was assessed over a 4-day time course after treatment with a non-targeting control shRNA (shCONT) or two independent non-overlapping shRNAs (shALKBH1.1008 and shALKBH1.1551) in two human patient-derived GSC models (387 and D456). Significance was determined by two-way (legend continued on next page) to non-tumor brain tissue and associated with reduced survival and increased glioma grade (Figures 7A-7C). An ALKBH1-regulated gene signature, defined by downregulated genes following ALKBH1 knockdown, correlated with tumor grade and shorter patient survival in several datasets ( Figures 7D-7G). Taken together, our findings demonstrate that ALKBH1 and regulation of N 6 -mA is necessary for glioma growth and may be a therapeutically targetable node.

DISCUSSION
N 6 -mA was originally described as a DNA modification that discriminates the original and newly synthesized DNA strand in bacteria (O'Brown and Greer, 2016). We and others reported that DNA undergoes rare N 6 -mA modification in low and high eukaryotes (Fu et al., 2015;Greer et al., 2015;Wu et al., 2016;Xiao et al., 2018;Yao et al., 2017;Zhang et al., 2015). However, the biological function of N 6 -mA in human cancers remains unclear. Here we elucidated the function of N 6 -mA in human glioblastoma. N 6 -mA levels were elevated in tumor relative to normal brain tissues and in GSCs compared with normal human astrocytes, suggesting that GSCs co-opt the regulatory potential of N 6 -mA to promote a tumorigenic chromatin landscape. As in murine ESCs, N 6 -mA repressed gene expression. Our ChIP-seq results strongly associated N 6 -mA sites with a common heterochromatin marker, H3K9me3, in glioblastoma, which has not been demonstrated previously. H3K9me3 is not involved in transcriptional silencing of N 6 -mA at young LINE-1 elements in ESCs. These results indicate that glioma cells have hijacked and further adapted N 6 -mA-mediated silencing mechanisms typically employed during early embryogenesis. Previous studies have implicated H3K9me3 in the regulation of DNA methylation and shown that complete removal of DNA methylation affects H3K9me3 levels and disrupts heterochromatin architecture (Du et al., 2015;Saksouk et al., 2014). Our findings corroborate the potential crosstalk between DNA N 6 -mA and H3K9me3 and possible collaboration of these modifications in transcription silencing in human cancers. Future experiments will inform the relationship between N 6 -mA, 5mC, and H3K9me3 to discern the temporal sequence of transcriptional silencing mediated by these marks. The precise roles of N 6 -mA and its interplay with the canonical 5mC DNA methylation mark as well as other histone modifications remain to be uncovered.
We demonstrated that ALKBH1 is a critical N 6 -mA demethylase that acts to regulate a number of gene networks critical to GSC identity, including the hypoxia response pathway. Most epigenetic regulation functions optimally in a specific range, with either loss or gain of function leading to disruption of cellular biology. Concordantly, elevation of N 6 -mA levels by targeting its demethylase ALKBH1 inhibited tumor formation. Although it may appear paradoxical, both N 6 -mA levels as well as the N 6 -mA demethylase were upregulated in glioblastoma. This apparent discrepancy can be resolved by understanding that the regulation of N 6 -mA is a highly dynamic process that depends on both a demethylase and a (currently unknown) methyltransferase. We hypothesize that the flux of N 6 -mA is important, not just the absolute levels. We interrogated the role of N 6 -mA in glioblastoma through modulating the only known regulator of N 6 -mA, ALKBH1. Disruption of the dynamic regulation of N 6 -mA through knockdown of ALKBH1 led to decreased proliferation, self-renewal, and tumor formation capacity, suggesting that efficient regulation of N 6 -mA is a GSC dependency.
Although we identified the N 6 -mA demethylase, further studies are required to identify the methyltransferases active in glioblastoma to model N 6 -mA in cancer. We posit that glioblastoma cells have elevated genome-wide N 6 -mA levels, whereas, at selected genomic loci, uncontrolled increases will be detrimental to the tumor cells. The active regulation of this modification at local sites by both a demethylase and methyltransferase likely plays an important role in the tuning of biological processes. Although global N 6 -mA enrichment is localized near genes involved in neurodevelopmental processes, modifications that are specifically regulated by ALKBH1 tend to be localized near hypoxia response genes. It should be noted that methyl-6-adenosine (m6A) modifications of RNA have been recently reported as well, including in glioblastoma (Cui et al., 2017;Fu et al., 2014;Zhang et al., 2017). The biochemical relationship between N 6 -mA and m6A modifications on DNA and RNA appear to be distinct.
In conclusion, we demonstrated that the N 6 -mA DNA modification exists in human tissues and is dynamically regulated in GSCs by the DNA demethylase ALKBH1 and that maintenance of this regulatory circuitry is critical for cell survival and proliferation. Our findings reveal a new class of DNA modification in human disease and identify a potential therapeutic vulnerability that can be exploited for cancer therapy. Therefore, our results repeated-measures ANOVA with Dunnett multiple test correction. p < 0.0001. Bottom: the knockdown efficiency of shRNAs targeting ALKBH1 was assessed by immunoblot. (B) Top: cell proliferation was assessed over a 4-day time course following treatment with a non-targeting control sgRNA (sgCONT) or two independent nonoverlapping sgRNAs (sgALKBH1#4 and sgALKBH1#5). Significance was determined by two-way repeated-measures ANOVA with Dunnett multiple test correction; p < 0.0001. Bottom: immunoblot showing the ALKBH1 protein level following treatment with a sgCONT and two independent non-overlapping sgRNAs targeting ALKBH1 (sgALKBH1#4 and sgALKBH1#5). (C) Tumor sphere formation efficiency and self-renewal capacity were measured by extreme in vitro limiting dilution assays (ELDAs) in two human patient-derived GSC models (387, left, and D456, right) after transduction with shCONT or shALKBH1. 387, p = 7.28e À14 ; D456, p = 4.19e À13 . (D) Kaplan-Meier curves depict survival of immunocompromised mice bearing intracranial tumors grown from human patient-derived GSC models (387 and D456) following transduction with shCONT or shALKBH1. Significance was determined by log rank analysis. **p < 0.01; n = 5 for each group. (E) Kaplan-Meier curve depicting survival of immunocompromised mice bearing intracranial tumors grown from human patient-derived GSC models (387) following transduction with sgRNA targeting ALKBH1 or a non-targeting control. Significance was determined by log rank analysis. **p < 0.01; n = 5 for each group. (A and B) Data are presented as mean ± SD. See also Figure S7.  (legend continued on next page) offer a novel therapeutic and discovery paradigm for deadly cancers.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
The authors declare no competing interests. Thorvaldsdó ttir, H., Robinson, J.T., and Mesirov, J.P. (2013). Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Brief. Bioinform. 14, 178-192. Turcan, S., Rohle, D., Goenka, A., Walsh, L.A., Fang, F., Yilmaz, E., Campos, C., Fabius, A.W., Lu, C., Ward, P.S., et al. (2012). IDH1 mutation is sufficient to establish the glioma hypermethylator phenotype. Nature 483, 479-483. were performed under an animal protocol approved by the University of California, San Diego Institutional Animal Care and Use Committee. Healthy, wild-type female mice of NSG (NOD.Cg-Prkdcscid Il2rgtm1Wjl/SzJ, The Jackson Laboratory, Bar Harbor, ME, USA) background, 4-6 weeks old, were randomly selected and used in this study for intracranial injection. Mice had not undergone prior treatment or procedures. Mice were maintained in 14 hours light/10 hours dark cycle by animal husbandry staff at the University of California, San Diego, with no more than 5 mice per cage. Housing conditions and animal status were supervised by a veterinarian. Animals were monitored until neurological signs were observed, at which point they were sacrificed. Neurological signs or signs of morbidity included hunched posture, gait changes, lethargy and weigh loss. In parallel survival experiments, mice were observed until the development of neurological signs. Healthy female mice of NSG, 4-6 weeks old, were randomly selected and used in this study for intracranial injection.

STAR+METHODS KEY RESOURCES
Additional experimental cell lines 293FT cells (ThermoFisher Scientific, Cat# R70007) were used to generate lentiviral particles as described in the method details section. 293FT cells were derived from embryonic kidney cells from a female human. Human astrocytes (ThermoFisher Scientific Cat#N7805100) were derived from an 18-week old female fetus and human astrocytes (SciencCell Cat#1800) were derived from a 23-week old human fetus of unknown gender (

METHOD DETAILS
Enzymatic hydrolysis of DNA The enzymatic digestion of DNA was performed as described here. One mg of genomic DNA was first spiked with uniformly 15 N-labeled 2 0 -deoxyadenosine and D 3 -labeled N 6 -methyl-2 0 -deoxyadenosine. DNA hydrolysis was performed with the addition of nuclease P1 (0.1 U), phosphodiesterase 2 (0.000125 U), erythro-9-(2-hydroxy-3-nonyl)adenine (EHNA) (0.25 nmol) and 4.5 mL solution containing 300 mM sodium acetate (pH 5.6) and 10 mM zinc chloride to the DNA solution. Incubation was carried out at 37 C for 24 hours. A second step digestion was carried out with 0.1 unit of alkaline phosphatase, 0.00025 unit phosphodiesterase I, 6 mL of 0.5 M Tris-HCl buffer (pH 8.9) at 37 C for 2 hours. Subsequently, the digestion mixture was neutralized by the addition of formic acid and the enzymes in the digestion mixture were removed by chloroform extraction.  Dot blotting DNA samples were denatured at 99 C for 10 minutes, cooled down on ice for 3 minutes, neutralized with 10% vol of 6.6 M ammonium acetate. Samples were spotted on the membrane (Amersham Hybond-N+, GE) and air dry for 5 minutes, followed by UV-crosslink (2 3 auto-crosslink, 1800 UV Stratalinker, STRATAGENE). Membranes were blocked in 5% nonfat dry milk in TBST for 2 hours at room temperature, incubated with N 6 -mA antibodies (1:1000, Synaptic Systems, 202-003) for 16 hours at 4 C. After 5 washes, membranes were incubated with HRP linked secondary anti-rabbit IgG antibody (1:5,000, Cell Signaling 7074) for 1 hour at room temperature. Signals were detected with ECL Plus Western Blotting Reagent Pack (GE Healthcare).

Immunofluorescent staining
Cells were fixed in 4% paraformaldehyde and permeabilized in PBS containing 0.3% Triton X-100 for 15 minutes. After permeabilization, samples were treated with 2N HCl for 30 min and subsequently neutralized for 10 minutes with 0.1 M sodium borate buffer pH 8.5. Then, the samples were blocked for 2h in blocking solution (5% donkey serum in PBS, RNaseA was added to a final concentration of 100 mg/ml), followed by incubation with N 6 -mA antibodies (1:200, Synaptic Systems, 202-003). Primary antibodies were incubated overnight at 4 C, followed by species appropriate secondary antibodies (1:1000, Alexa 488; Invitrogen Molecular Probes, Eugene, OR) with incubation for 1 hour. Nuclei were stained with DAPI, and slides were then mounted using Fluoromount (Calbiochem, San Diego, CA). Images were taken using a Leica DM4000 Upright microscopy.

Proliferation assay
Cell proliferation was measured using Cell-Titer Glo (Promega, Madison, WI). For this assay, 1,000 cells were seeded into each well of a 96-well plate. Cells were plated in 100 mL of Neurobasal media (Invitrogen) supplemented with B27 without vitamin A (Invitrogen), EGF, and bFGF (20 ng/ml each; R&D Systems), sodium pyruvate, and glutamax. To assess cell viability, 40 mL of the Cell-Titer Glo (Promega, Madison, WI) reagent was added to each well. Plates were placed on an orbital shaker at 100rpm for 10 minutes to promote cell lysis and covered with aluminum foil to prevent exposure to light. Plates were incubated in the dark for 2 minutes and then plates were read using a luminometer. All data were normalized to day 0 and presented as mean ± SEM.
Generation of lentiviral particles 293FT cells (ThermoFisher Scientific, Cat# R70007) were used to generate lentiviral particles through co-transfection of the packaging vectors pCMV-dR8.2 dvpr (Addgene, Plasmid # 8455) and pCl-VSVG (Addgene, Plasmid # 1733) with the appropriate overexpression, shRNA, or other plasmid. 293FT cells were seeded at a density of 1.2 million cells in DMEM, high glucose (ThermoFisher Scientific, Cat# 11995073) in 10% Fetal Bovine Serum (ThermoFisher Scientific, Cat# 26140079 with 1% Penicillin-Streptomycin (ThermoFisher Scientific, Cat # 15140122). Cells were incubated for 24 hours prior to transfection. Transfection was performed using LipoD293 In Vitro DNA Transfection Reagent (SignaGen Laboratories, Cat # SL100668) according to the manufacturer's instructions. Briefly, 5 mg each of the packaging plasmids and the plasmid of interest were combined into a tube, the transfection reagent was diluted and added followed by a 15 minute incubation. The transfection mixture was then added to the 293FT cells. Media was changed after 16 hours. Virus was collected 48 hours after media change and concentrated using the Lenti-X Concentrator (Clontech Takara Bio USA, Cat # 631232) according to the manufacturers instructions. Viral supernatants were centrifuged at 1,500 x g for 45 minutes and viral pellets were frozen at À80 C for future use.

In vitro limiting dilution assay
For in vitro limiting dilution assays, decreasing numbers of cells per well (20, 10, 5, and 1) were plated in 96-well plates in neurobasal media supplemented with the components listed above. Ten days after plating, the presence and number of neurospheres in each well was quantified. Extreme limiting dilution analysis was performed using software available at http://bioinf.wehi.edu.au/software/ elda, as previously described.

Glioblastoma Stem Cell Differentiation
For in vitro differentiation, glioblastoma stem cells were cultured for one week in DMEM, high glucose (ThermoFisher Scientific, Cat# 11995073) in 10% Fetal Bovine Serum (ThermoFisher Scientific, Cat# 26140079 with 1% Penicillin-Streptomycin (ThermoFisher Scientific, Cat # 15140122) as part of a previously established differentiation protocol. Cellular differentiation was verified by observing depleted expression of SOX2 protein by western blot.

Intracranial tumor formation and in vivo bioluminescence imaging
GSCs were transduced with lentiviral vectors expressing ALKBH1 or a non-targeting, control (shCONT) shNRA for the knockdown experiments. 36 hours post infection, viable cells were counted and engrafted intracranially into NSG (NOD.Cg-Prkdcscid Il2rgtm1Wjl/SzJ, The Jackson Laboratory, Bar Harbor, ME) mice under a University of California, San Diego Institutional Animal Care and Use Committee (IACUC) approved protocol. In parallel survival experiments, animals were monitored until they developed neurological signs.
N 6 -mA DIP-seq Genomic DNA from patient-derived glioblastoma stem cell models was purified with a DNeasy kit (QIAGEN,69504). For each sample, 5 mg DNA was sonicated to 200-300 bp with Bioruptor. Then, adaptors were ligated to genomic DNA fragments following the Illumina protocol. The ligated DNA fragments were denatured at 95 degree for 5 minutes. Then, the single-stranded DNA fragments were immunoprecipitated with N 6 -mA antibodies (5 mg for each reaction, 202-003, Synaptic Systems) overnight at 4 C. N 6 -mA enriched DNA fragments were purified according to the Active Motif hMeDIP protocol. IP DNA and input DNA were PCR amplified with NEBNext indexing primers, and were then subjected to multiplexed library construction and sequencing with Illumina HiSeq sequencing.

ChIP-sequencing
Formaldehyde-fixed cells were lysed and sheared (Branson S220) on wet ice. The sheared chromatin was cleared and incubated overnight at 4 C with 5 mg H3K9me3 antibody (Abcam, ab8898), H3K27me3 (active motif, 39155), or H3K4me3 (Abcam) respectively. Antibody-chromatin complexes were immunoprecipitated with protein G magnetic Dynal beads (Life Technologies), washed, eluted, reverse crosslinked, and treated with RNase A followed by proteinase K. ChIP DNA was purified using Ampure XP beads (Beckmann Coulter) and then used to prepare sequencing libraries by NEB NEBNextâ Ultra DNA Library Prep Kit for Illumina for sequencing.

ChIP-sequencing Data Analysis
Raw FASTQ files were trimmed using TrimGalore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) and aligned to the hg19 human genome using the Burrows-Wheeler Aligner (BWA) (Zhang et al., 2008). Identical ChIP-seq sequence reads were collapsed to a single read to avoid PCR duplicates using Picard and SamTools was used to sort and index the BAM files. N 6 -mA Specific Analysis Parameters BAM files from the N 6 -mA DIP-seq and matched inputs were converted to BED using the Bedtools ''bamtobed''function (Quinlan and Hall, 2010). N 6 -mA peaks were called with the SICER algorithm (Xu et al., 2014) using the EPIC wrapper script (https://github.com/ biocore-ntnu/epic) with input DIP-seq files used as a control and using the default significance cutoff of False Discovery Rate (FDR) < 0.05. Peak intersections were calculated using the Bedtools ''intersect'' function (Quinlan and Hall, 2010). HOMER was used to perform genome annotation using the ''annotatePeaks'' function and for motif analysis using the ''findMotifsGenome'' function (Heinz et al., 2010). CEAS was also used for genome ontology analysis and for determining chromosome distribution of DIPseq peaks (Ji et al., 2006). For gene set enrichment analysis, peaks were mapped to the closest transcription start site using the HOMER ''annotatePeaks'' function. Genes were then used as inputs into the online Gene Set Enrichment Analysis web portal (http://software.broadinstitute.org/gsea/index.jsp) (Subramanian et al., 2005) (Mootha et al., 2003. Pathway enrichment bubble plots were generated using the Bader Lab Enrichment Map Application (Merico et al., 2010) and Cytoscape (https://cytoscape.org). For subsequent visualization of N 6 -mA signal, bigwig files were generated using the DeepTools ''bamCoverage'' script using default settings (https://deeptools.readthedocs.io/) (Ramírez et al., 2014). Signal tracks were visualized using the Integrative Genomics Viewer, IGV (Robinson et al., 2011;Thorvaldsdó ttir et al., 2013). For the generation of heatmaps and profile plots, the DeepTools ''compu-teMatrix,'' ''plotHeatmap,'' and ''plotProfile'' functions were used with specific parameters given in the figure legends of the manuscript.
Histone Modification Specific Analysis Parameters BAM files from H3K9me3, H3K27me3, and H3K4me3 ChIP-seq experiments and matched inputs were converted to BED using the Bedtools ''bamtobed '' function (Quinlan and Hall, 2010). Histone modification peaks were called with the SICER algorithm (Xu et al., 2014) using the EPIC wrapper script (https://github.com/biocore-ntnu/epic) with input ChIP-seq files used as a control and using the default significance cutoff of False Discovery Rate (FDR) < 0.05. Intersections between N 6 -mA and various histone modification peaks were calculated using the Bedtools ''intersect'' function (Quinlan and Hall, 2010). To identify regions of gained H3K9me3 enrichments following ALKBH1 knockdown, the SICER algorithm was used as described above with the ALKBH1-knockdown H3K9me3 ChIP-seq file serving as the ''treatment'' file and the non-targeting shRNA H3K9me3 ChIP-seq file serving as the ''control'' file using default settings as described above. Reciprocally, to identify regions of lost H3K9me3 enrichments following ALKBH1 knockdown, the SICER algorithm was used as described above with the non-targeting H3K9me3 ChIP-seq file serving as the ''treatment'' file and the ALKBH1-knockdown H3K9me3 ChIP-seq file serving as the ''control'' file. Consensus gained and lost peaks were generated by identifying peaks observed in both technical replicates using the bedtools ''intersect'' function (Quinlan and Hall, 2010) and results are plotted as a volcano plot using the R programming language. For subsequent visualization of histone modification ChIP-seq signal, bigwig files were generated using the DeepTools ''bamCoverage'' script using default settings (https:// deeptools.readthedocs.io/) (Ramírez et al., 2014). Signal tracks were visualized using the Integrative Genomics Viewer, IGV (Robinson et al., 2011;Thorvaldsdó ttir et al., 2013). For the generation of heatmaps and profile plots, the DeepTools ''computeMatrix,'' ''plotHeatmap,'' and ''plotProfile'' functions were used with specific parameters given in the figure legends of the manuscript. ALKBH1 ChIP-seq Specific Parameters Peaks were called using the MACS2 ''callpeak'' function using default settings using paired-end BAM files from ALKBH1 ChIP-seq experiments using the ''-f BAMPE'' parameter and the default FDR-corrected q-value of 0.05 with a matched ChIP-seq input file serving as the control. The bedtools ''intersect'' function was used to calculate overlaps with N 6 -mA peaks, and the bedtools ''shuffle'' function was used to generate a ''genome background'' expected overlaps file (Quinlan and Hall, 2010). For subsequent visualization of ALKBH1 ChIP-seq signal, bigwig files were generated using the DeepTools ''bamCoverage'' script using default settings (https://deeptools.readthedocs.io/) (Ramírez et al., 2014). For the generation of heatmaps and profile plots, the DeepTools ''computeMatrix,'' ''plotHeatmap,'' and ''plotProfile'' functions were used with specific parameters given in the figure legends of the manuscript.

RNA-sequencing
Cellular RNA was extracted using the QIAGEN miRNeasy Mini Kit (catalog number 217004) according to the manufacturers instructions. Total RNA was prepared for sequencing using the Illumina TruSeq Stranded Total RNA Library Prep Kit. For samples treated in normoxia, FASTQ sequencing files were trimmed and then aligned to the hg19 human genome using the STAR aligner (Dobin et al., 2013). Differential expression analysis was performed using CuffDiff (http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/). For samples treated in hypoxia, FASTQ sequencing reads were trimmed using Trim Galore (https://www.bioinformatics.babraham.ac.uk/ projects/trim_galore/) and transcript quantification was performed using Salmon in the quasi-mapping mode (Patro et al., 2017). Salmon ''quant'' files were converted using Tximport (https://bioconductor.org/packages/release/bioc/html/tximport.html) and differential expression analysis was performed using DESeq2 (Love et al., 2014). Gene set enrichment analysis was performed by selecting differentially expressed genes (FDR-corrected p value < 0.05), generating a pre-ranked list using the gene expression fold change as the ranking metric, and inputting a pre-ranked list into the GSEA desktop application (

Bisulfite Sequencing
Genomic DNA was extracted from human patient-derived glioblastoma stem cells and bisulfite conversion was performed. Following paired-end sequencing, FASTQ files separated into 11 paired-end fastq files to facilitate processing. Reads were trimmed using Trim Galore to remove 10bp from the 5 0 and 3 0 end of both reads. Bowtie2 (Langmead and Salzberg, 2012) was used to map the reads to a pre-prepared hg19 genome in the non-directional and paired-end mode. The De-duplicate function in Bismark (Krueger and Andrews, 2011) was used to remove duplicate reads. All of the mapping results from each pair of small fastq files were merged together using the samtools merge method, with the -n option to sort the reads by names. The function ''bismark_methylation_extractor'' was applied on the merged mapping results to extract the methylation information. The coverage2cytosine function was used to generate a genome-wide cytosine methylation report. The MakeTagDirectory function in homer was applied on the cytosine methylation report, with the option of -format bismark, -minCounts 5 -checkGC to calculate percentage of methylation in the CpG locations across the genome. The annotatePeaks.pl script from HOMER (Heinz et al., 2010) was used to obtain the methylation level of selected regions of interest. We calculated the average methylation in the 12,000 bp regions centered on the N 6 -mA peaks as described in the manuscript and TSS. The cytosine methylation report was also used to calculate the methylation level for each individual N 6 -mA peak, 1000 bp regions ahead of the TSS, and also for a bed file with shuffled chip-seq peak across the genome. Three violin plots were generated. Rank-sum test was applied between three of the distributions.

QUANTIFICATION AND STATISTICAL ANALYSIS
Analysis for each plot is listed in the figure legend and/or in corresponding methods above. In brief, all grouped data are presented as mean ± s.d. All box and whisker plots of expression data are presented as median (middle line of box) ± the 25 percentile (top and bottom line of box, respectively). P values presented are calculated by two-sided Mann-Whitney U-test. Kaplan-Meier curves were generated and log-rank (Mantel-Cox) analysis was performed to generate P values using GraphPadPrism software (GraphPad Software). Sample sizes for each experiment are given in corresponding figures and/or methods above. Sizes were chosen based on previous experience with given experiments, or in the case of retrospective analysis, all available samples were included.

DATA AND SOFTWARE AVAILABILITY
All raw and selected processed data files are available on the Gene Expression Omnibus (GEO). All data can be accessed at the SuperSeries accession number GEO: GSE118093 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE118093). Direct links to the RNA-sequencing data can be found at the SubSeries accession number GEO: GSE117632. Direct links to ATAC-sequencing data can be found at the SubSeries accession number GEO: GSE118092. Direct links to ChIP-sequencing, whole genome bisulfite sequencing, and N 6 -mA DIP-sequencing data can be found at the SubSeries accession number GEO: GSC119081. Additional data will be provided upon request. There are no restrictions on data availability.      A) Visualization of N 6 -mA DIP-seq and H3K9me3 ChIP-seq signal surrounding the CDKN3 locus. (B) Visualization of N 6 -mA DIP-seq and H3K9me3 ChIP-seq signal surrounding the RASSF2 locus. (C) Patients in the TCGA HG-U133A microarray database were stratified based on median expression of AKAP6 into ''high'' (greater than or equal to the median expression, n = 263) or ''low'' (less than median expression, n = 262) groups. Survival is shown as a Kaplan-Meier plot. Significance was determined by log-rank test. p = 0.0316. (D) Log mRNA expression of AKAP6 in non-tumor tissue and glioblastoma tissue is plotted from the TCGA HG-U133A microarray database. Significance was determined by an unpaired two-tailed t test. p < 0.001. (E) Visualization of N 6 -mA DIP-seq and H3K9me3 ChIP-seq signal surrounding the AKAP6 locus.