Abstract
Glioblastoma is a universally lethal cancer driven by glioblastoma stem cells (GSC). Here, we interrogated N6-methyladenosine (m6A) mRNA modifications in GSCs by methyl RNA immunoprecipitation followed by sequencing and transcriptome analysis, finding transcripts marked by m6A often upregulated compared with normal neural stem cells (NSC). Interrogating m6A regulators, GSCs displayed preferential expression, as well as in vitro and in vivo dependency, of the m6A reader YTHDF2, in contrast to NSCs. Although YTHDF2 has been reported to destabilize mRNAs, YTHDF2 stabilized MYC and VEGFA transcripts in GSCs in an m6A-dependent manner. We identified IGFBP3 as a downstream effector of the YTHDF2–MYC axis in GSCs. The IGF1/IGF1R inhibitor linsitinib preferentially targeted YTHDF2-expressing cells, inhibiting GSC viability without affecting NSCs and impairing in vivo glioblastoma growth. Thus, YTHDF2 links RNA epitranscriptomic modifications and GSC growth, laying the foundation for the YTHDF2–MYC–IGFBP3 axis as a specific and novel therapeutic target in glioblastoma.
Epitranscriptomics promotes cellular heterogeneity in cancer. RNA m6A landscapes of cancer and NSCs identified cell type–specific dependencies and therapeutic vulnerabilities. The m6A reader YTHDF2 stabilized MYC mRNA specifically in cancer stem cells. Given the challenge of targeting MYC, YTHDF2 presents a therapeutic target to perturb MYC signaling in glioblastoma.
This article is highlighted in the In This Issue feature, p. 211
Introduction
Glioblastoma (World Health Organization grade IV glioma) represents the most common primary, intrinsic brain tumor with inevitable recurrence, limiting the median survival of patients to little more than a year (1, 2). Glioblastomas display cellular hierarchies with self-renewing glioblastoma stem cells (GSC) at the apex, with contributions of GSCs to therapeutic resistance and tumor recurrence (3–5). Standard-of-care therapy includes surgical resection followed by combined radiotherapy and chemotherapy, and then adjuvant chemotherapy, but treatment remains palliative (6). Given the roles of GSCs in therapeutic resistance, angiogenesis, immune escape, and invasion, clinical and preclinical observations suggest that targeting GSCs may improve tumor outcome (7).
Precision medicine efforts in neuro-oncology have demonstrated limited benefit in most patients, prompting interrogation of epigenetic cell states to inform therapeutic targeting (8). Posttranscriptional regulation of RNA, epitranscriptomics, adds to the layers of regulation in normal and neoplastic cell biology. N6-methyladenosine (m6A) is the most prevalent eukaryotic RNA modification (9) and regulates the stability and translation of modified mRNAs (10–13). Like other methylation events, m6A is dynamically regulated by three classes of proteins: a “writer complex” that installs m6A marks (methyltransferase components: METTL3, METTL14, and WTAP); “erasers” that remove m6A marks (demethylases: FTO and ALKBH5); and “readers” and other RNA-binding proteins that recognize and bind to m6A-modified transcripts and mediate the downstream signaling events (m6A-binding proteins: YTHDF1–3 and YTHDC1–2; refs. 14, 15). m6A writers (16–18) and m6A erasers (19–21) are differentially regulated in cancers, including glioblastoma, but the downstream connection between RNA m6A methylation and tumor biology is poorly understood. Recently, m6A readers were implicated in tumor biology by contributing to tumor growth, proliferation, and survival (22–24). m6A reader proteins contain a YTH domain (YT521-B homology), which mediates specific RNA binding (25). There are five YTH domain–containing proteins, including YTHDC1, which is predominantly nuclear and regulates posttranscriptional splicing events (26), and YTHDC2 and YTHDF1–3, which are predominantly cytosolic. Some studies have suggested that recognition of m6A by YTHDF1 and YTHDF3 facilitates protein translation (27, 28), whereas YTHDF2 and YTHDF3 transport mRNAs from the actively translating pool of mRNAs to the decaying mRNA pool, leading to mRNA degradation (29). In contrast, more recent reports suggest that all three YTHDF proteins function to mediate degradation of m6A-tagged mRNAs (13, 30). Here, we compared the transcriptomic m6A landscapes of normal and neoplastic stem cells, investigated the roles of m6A readers in GSCs and glioblastoma, and interrogated the mechanism through which these epitranscriptomic modifications support tumor biology. We report that YTHDF2 stabilizes important oncogenic drivers in an m6A-dependent manner specific to GSCs.
Results
Oncogenic Transcripts Upregulated in GSCs Are Marked by RNA m6A Modification
RNA methylation is one of the most prevalent RNA modifications in human cells (31). m6A decorates messenger RNAs, long noncoding RNAs, and ribosomal and spliceosome RNAs (32). To understand the global distribution of m6A marks, we performed transcriptome-wide m6A methylation analyses in two patient-derived GSCs (4121 and 387) and two neural stem cells (NSC; NSC11 and HNP1; Supplementary Fig. S1A). m6A signals were normalized to respective inputs to adjust for baseline transcript expression levels. m6A peaks discriminated GSCs from NSCs in unsupervised clustering analyses (Supplementary Fig. S1B), suggesting that GSCs demonstrate altered m6A distribution compared with their nonneoplastic counterparts, confirmed by the identification of specific m6A peaks gained in GSCs relative to NSCs (Fig. 1A; Supplementary Fig. S1C; Supplementary Table S1). Motif analysis using the top 1,000 enriched peaks from GSCs and NSCs identified highly enriched, conserved consensus motifs in both cell types (Fig. 1B). m6A peak density was not globally altered between NSCs and GSCs, with approximately equal numbers of increased and decreased m6A peaks in each cell type (Supplementary Fig. S1D). Total m6A levels were quantified by ELISA after extracting total RNA (which reflects m6A from rRNA, small nuclear RNA, and mRNA) from six functionally validated human patient–derived GSCs (387, 4121, GSC23, 3565, 738, and 1919) and four human NSCs (NSC11, HNP1, ENSA, and NSC194). Total m6A RNA levels did not significantly differ between GSCs and NSCs (Supplementary Fig. S1E). m6A peaks were distributed throughout transcripts, with the majority found on the gene body (Supplementary Fig. S1F). Analysis of m6A signal distribution along mRNA transcripts confirmed enrichment of m6A signals around stop codons, as previously reported (33, 34), along with the strongest enrichment near start codons in all cell types (Supplementary Fig. S1G).
We next investigated whether differential localization of m6A peaks informed GSC biology. We generated a GSC-specific gene expression signature by comparing RNA sequencing (RNA-seq) data from our cohort of 38 GSCs and 5 NSCs (8). Gene set enrichment analysis (GSEA) showed that genes with gained m6A peaks were highly enriched with the GSC signature (Fig. 1C). To assess the effect of m6A marks on gene expression, we leveraged the consensus m6A peak set to identify genes with or without m6A marks. Transcripts with m6A peaks displayed higher expression than transcripts lacking m6A marks (Supplementary Fig. S1H and S1I). This may in part reflect improved detection of m6A on more highly expressed mRNAs (35). We found genes with m6A peaks gained in GSCs were globally upregulated (Fig. 1D and E). Conversely, genes with m6A peaks lost in GSCs relative to NSCs were generally downregulated in GSCs (Fig. 1F). As detection of m6A peaks on genes expressed at low levels is technically challenging, we filtered out genes expressed at low levels across the RNA-seq cohort to minimize bias. In GSCs, m6A peaks were gained on important genes associated with cancer stem cells, including OLIG2 and MYC, concurrent with increased expression (Fig. 1G and H). MN1 [meningioma (disrupted in balanced translocation) 1] showed the highest gain in m6A signal in NSCs and was expressed at much higher levels in NSCs compared with GSCs (Fig. 1I and J). Although these changes in m6A may be due to changes in the underlying mRNA expression, these changes may also reflect a stabilizing role for m6A on mRNA in these cells.
To validate these observations in an independent dataset, we mined the pan-glioma dataset of The Cancer Genome Atlas (TCGA; ref. 36) to explore mRNA expression levels of m6A-marked genes in patient tumors. MYC, which gained multiple m6A peaks in GSCs, was overexpressed in gliomas (Fig. 1K), whereas MN1, which lost m6A peaks in GSCs, was downregulated in gliomas (Fig. 1L). The top 100 genes marked with the highest m6A signals in GSCs were overexpressed in TCGA glioblastoma tumors versus normal brain (Fig. 1M), suggesting that in GSCs, m6A modification marks important genes overexpressed in glioblastomas. Genes marked by GSC-enriched m6A peaks were ranked based on both statistical significance predicting survival of patients afflicted with glioblastoma and differential duration in survival when overexpressed, relative to median mRNA levels. Several genes, including NOL3, EMILIN1, and TRAF1, portend poor prognosis when overexpressed in tumors of patients with glioblastoma (Supplementary Fig. S1J), suggesting that m6A peaks mark potential glioblastoma oncogenes. Furthermore, GSC-enriched m6A marks were present on transcripts essential for GSC viability by whole-genome CRISPR screening (37), including MYC, ATP2A2, and IMP4 (Supplementary Fig. S1K). Subsetting transcripts by the localization of m6A peaks revealed that the location of m6A marks was associated with differential transcript expression. Transcripts with m6A peaks on the 5′ untranslated region (UTR) or promoter region (which may also reflect m6Am rather than m6A) were expressed at higher levels than those harboring peaks on the 3′ UTR, gene body, or stop codon (Supplementary Fig. S1L and S1M). Collectively, these results indicate that both the presence and localization of m6A may affect gene expression. We performed pathway enrichment analysis on gene sets marked with cell type–specific m6A peaks. Genes with GSC-enriched m6A peaks were enriched for programs of cell communication and signal transduction, transcriptional regulation, and downregulation of UV response, among others (Fig. 1N), whereas genes with NSC-enriched m6A peaks were selectively enriched for developmental pathways, including tissue and nervous system development (Fig. 1O). Taken together, these results implicate m6A peaks in marking, and potentially regulating, critical cancer-specific genes in glioblastoma in a location and cell type–specific manner.
YTHDF2 Is Upregulated in GSCs and Essential for GSC Maintenance
As m6A deposition and localization may regulate cell type–specific processes in GSCs, we investigated the relative contributions of m6A writers, erasers, and readers in glioblastoma. In TCGA glioma patient tumor data, expression of m6A writers (METTL3 and METTL14) or m6A erasers (ALKBH5 and FTO) did not increase with tumor grade (Fig. 2A; Supplementary Fig. S2A and S2B). Of the m6A readers, expression of YTHDF2 and YTHDF3 was upregulated in high-grade gliomas (Fig. 2A). Further, YTHDF2 expression was enriched in mesenchymal glioblastomas (Supplementary Fig. S2B), which demonstrate greater therapeutic resistance and worse prognosis (38, 39). We compared expression of m6A regulators in GSCs versus NSCs; GSCs expressed higher levels of the m6A writer METTL14 and the m6A reader YTHDF2 than NSCs, whereas other m6A writers and erasers were not differentially expressed between these cell types (Fig. 2B and C). To investigate the expression of METTL14 and YTHDF family members in the absence of cell culture conditions, we mined the pan-glioma TCGA dataset. Expressions of YTHDF family members (YTHDF1–3; Supplementary Fig. S2C) and METTL14 (Supplementary Fig. S2D) were elevated in glioblastoma tissues compared with the normal brain tissue.
As YTHDF2 and METTL14 expression distinguishes GSCs from NSCs and glioblastoma from normal brain, we investigated the functional roles of these molecules both in vitro and in vivo. METTL14 promotes leukemogenesis through the maintenance of hematopoietic stem cells (18). Targeting METTL14 with CRISPR/Cas9 impaired the proliferation capacity of GSCs in vitro (Supplementary Fig. S2E and S2F), but did not affect the growth of intracranial glioblastoma tumors derived from human patient–derived GSCs (Supplementary Fig. S2G), suggesting a limited role in cancer stemness.
We then investigated the functional role of the m6A YTHDF readers in glioblastoma. We first examined YTHDF1 and YTHDF3, which are typically expressed at much lower levels than YTHDF2 in most cells (40). Targeting either YTHDF1 or YTHDF3 expression with two independent single guide RNAs (sgRNA) per target gene had minimal effect on the proliferation of patient-derived GSCs relative to a control nontargeting sgRNA (Fig. 2D–F). In contrast, knockout of YTHDF2 with either of two independent sgRNAs decreased cell viability relative to a control sgRNA (Fig. 2G and H). GSC frequency and self-renewal were diminished, as assayed by extreme limiting dilution assays (ELDA), following depletion of YTHDF2 in two different GSC models (Fig. 2I). Sphere formation under serum-free conditions has served as a surrogate for self-renewal, albeit with caveats (41). YTHDF2 knockout reduced sphere formation in two GSCs relative to the control sgRNA (Fig. 2J). We next investigated whether YTHDF2 depletion induces GSC differentiation. Upon YTHDF2 knockout in GSCs, expression of stemness markers SOX2 and OLIG2 decreased without change in differentiation marker GFAP (Supplementary Fig. S2H and S2I). YTHDF2 functions were not selective to GSCs alone, as YTHDF2 depletion had similar effects on cell viability in matched GSCs and differentiated glioma cells (Supplementary Fig. S2J and S2K). To verify these phenotypes, we performed orthogonal studies using two independent shRNAs targeting YTHDF2 with rescue overexpression of an shRNA-resistant YTHDF2 (Fig. 2K). Concordant with the CRISPR studies, shRNA-mediated knockdown of YTHDF2 reduced GSC viability, which was rescued upon overexpression of shRNA-resistant YTHDF2 (Fig. 2L). YTHDF2 overexpression in GSCs was confirmed by using immunofluorescence analysis (Fig. 2M), which revealed predominantly cytosolic localization of YTHDF2. Collectively, these results establish YTHDF2 as a specific and potent regulator of glioblastoma maintenance.
YTHDF2 Supports Gene Expression in GSCs through m6A RNA Modification
To determine the mechanism by which YTHDF2 functions in GSCs, we interrogated the downstream targets of YTHDF2 using RNA-seq. Knockout of YTHDF2 induced widespread gene-expression changes in GSCs (Fig. 3A; Supplementary Table S2). Programs enriched upon YTHDF2 knockdown included translational control, ribosomes and RNA processing, hypoxia response, and MYC targets (Fig. 3B and C). Transcriptional changes following YTHDF2 knockdown aligned with differential m6A peaks between GSCs and NSCs. Genes with gained m6A peaks in GSCs were more frequently downregulated upon YTHDF2 knockout, supporting YTHDF2 regulation of GSC-specific transcriptional programs in an m6A-dependent manner (Fig. 3D). Targeting YTHDF2 reduced expression of MYC targets and genes involved in hypoxia, ribosomes, glycolysis, RNA processing, and translation (Supplementary Fig. S3A). The effects of YTHDF2 knockout on decreasing mRNA levels of MYC, VEGFA, and other targets (including IGFBP3, IGFBP5, DDIT, and ribosomal genes RPS2 and RPL13A) were validated by qPCR (Fig. 3E and F). Targeting YTHDF2 reduced protein levels of the GSC stemness markers SOX2 and OLIG2 (Supplementary Fig. S3B), connecting YTHDF2 depletion and loss of GSC maintenance (Fig. 2I). MYC (42) and IGFBP3 (43) have been implicated in maintenance of tissue and cancer stem cell function. The lower transcript levels of MYC and IGFBP3 upon YTHDF2 knockout translated into decreased protein levels in GSCs, measured by immunoblot (Supplementary Fig. S3B). To predict the function of YTHDF2 in glioblastoma, we interrogated TCGA glioblastoma gene-expression data and identified programs correlated with YTHDF2 expression. YTHDF2-correlated genes were positively enriched for processes of RNA metabolism and translation, cell cycle, glucose metabolism, and DNA repair signatures, whereas YTHDF2 expression negatively correlated with IFN and GPCR signaling (Fig. 3G). YTHDF2-correlated genes were highly enriched for MYC and E2F targets and for G2–M checkpoint regulators and mediators of oxidative phosphorylation (Fig. 3H). These data implicate YTHDF2 as a regulator of global transcriptional programs associated with differential m6A modification.
To consider YTHDF2 and MYC in glioblastoma subtypes, we explored the correlation between YTHDF2 and MYC mRNA expression in all glioblastoma specimens in an independent dataset, the Chinese Glioma Genome Atlas (CGGA). YTHDF2 mRNA expression correlated with MYC across glioblastoma in general, as well as proneural, mesenchymal, and classic subtypes individually (Supplementary Fig. S3C–S3F), suggesting that the correlation between MYC and YTHDF2 in our patient-derived GSCs is consistent across glioblastoma transcriptional subtypes.
YTHDF2 Stabilizes MYC Transcripts
Following YTHDF2 knockout, mRNA transcripts were upregulated (n = 498) and downregulated (n = 512) in almost equal proportion (Fig. 3A). As YTHDF2 is believed to bind and promote degradation of m6A-tagged transcripts (29), we asked whether genes perturbed upon YTHDF2 knockout included direct targets of YTHDF2. To identify the direct targets of YTHDF2, we mapped the binding sites of YTHDF2 using individual nucleotide resolution cross-linking and immunoprecipitation (iCLIP) in two patient-derived GSCs (4121 and 387) and two NSCs (NSC11 and HNP1). YTHDF2 binding to RNAs in NSC11 was very low, precluding library preparation and sequencing in this NSC line (Supplementary Fig. S4A), likely as a result of low YTHDF2 levels in these cells. Motif analysis using the top enriched peaks from the two GSCs and HNP1 NSCs identified highly enriched, conserved consensus motifs in both cell types (Supplementary Fig. S4B). Motifs identified were similar to m6A consensus DRACH motifs (Supplementary Fig. S4B), indicating that YTHDF2 bound to m6A peak regions in these cells. YTHDF2-binding sites were distributed throughout transcripts with an enrichment around the stop codon and 3′ UTR region in each cell type (Supplementary Fig. S4C).
We hypothesized that critical targets of YTHDF2 would both bind to YTHDF2 (as measured by YTHDF2 iCLIP-seq) and be reduced in expression upon YTHDF2 targeting. Upon overlap of YTHDF2 iCLIP-seq data and YTHDF2-knockout RNA-seq data from GSCs, we identified 63 genes directly bound by YTHDF2 among 512 genes downregulated upon YTHDF2 knockout (Fig. 4A; Supplementary Table S3). This restricted gene set enriched for pathways associated with translation initiation, ribosome complexes, nonsense-mediated decay, and carbon metabolism in cancer (Fig. 4B). To further refine the downstream targets, we focused on the overlap between transcripts that bound YTHDF2 and had detectable m6A marks among genes downregulated upon YTHDF2 targeting, identifying 30 genes, including important oncogenes such as MYC (42, 44) and VEGF (ref. 45; Fig. 4C). In contrast, YTHDF2 bound to 33 transcripts that were upregulated upon YTHDF2 knockout (Supplementary Fig. S4D; Supplementary Table S3), of which 23 were m6A-marked (Supplementary Fig. S4E).
YTHDF2 displayed selective binding to MYC mRNA in patient-derived GSCs, but not in HNP1 NSCs (Fig. 4D). To validate our findings in GSCs, we analyzed available YTHDF2 photoactivatable ribonucleoside–enhanced CLIP data and m6A profiling data (29) obtained from HeLa cells and overlapped these results with our YTHDF2-knockout RNA-seq data obtained from GSCs. In agreement with our results, this study reported the presence of m6A and direct binding of YTHDF2 to MYC and VEGFA in HeLa cells (Supplementary Fig. S5A–S5C). YTHDF2-bound genes in HeLa cells were enriched for pathways associated with protein translation initiation and glycolysis, providing an independent validation to our findings (Supplementary Fig. S5A and S5C).
We validated YTHDF2 binding to its targets in GSCs by YTHDF2 RNA immunoprecipitation (RIP) coupled with qPCR. YTHDF2 binding was selectively enriched on MYC and VEGF transcripts in GSCs compared with NSCs (Fig. 4E and F). YTHDF2 binding to CBP mRNA was used as a positive control (Supplementary Fig. S5D). To determine whether YTHDF2 regulates the stability of these transcripts, we treated GSCs transduced with either control sgRNA or YTHDF2 sgRNA with actinomycin D to arrest transcription. The decay rates of MYC and VEGFA mRNAs were higher upon YTHDF2 depletion, suggesting that YTHDF2 is critical for the stabilization of these transcripts (Fig. 4G and H). Taken together, our data raise the possibility that YTHDF2 binding may not only destabilize certain transcripts (24, 29), as shown by mRNAs that were m6A-tagged and upregulated upon YTHDF2 knockdown (Supplementary Fig. S4D and S4E), but also stabilize important oncogenes, such as MYC and VEGF, in an m6A-dependent manner.
To further establish the role of m6A in YTHDF2 functions, we knocked out the m6A writer METTL3 using CRISPR/Cas9 (Supplementary Fig. S5E) and performed RIP-qPCR to measure YTHDF2 binding to MYC mRNA. As expected, METTL3 knockout decreased total m6A levels in GSCs (Fig. 4I). METTL3 knockout reduced YTHDF2 binding to MYC mRNA (Fig. 4J and K) and decreased MYC mRNA levels in GSCs (Supplementary Fig. S5F). To address the importance of m6A for YTHDF2-mediated transcript stabilization, we generated FLAG-tagged m6A-binding mutant YTHDF2 constructs (W432A and W486A) and overexpressed wild-type and m6A-binding mutant YTHDF2 in GSCs. Using RIP-qPCR, we found reduced binding of mutant YTHDF2 to MYC mRNA in comparison with wild-type YTHDF2 (Fig. 4L and M). These results confirmed that m6A was necessary for the effects of YTHDF2 on MYC expression. To confirm that the m6A-binding function of YTHDF2 was critical to its role as a GSC dependency, we rescued YTHDF2 knockdown cells with either wild-type or m6A-binding mutant YTHDF2 constructs. Only wild-type YTHDF2 rescued the viability of YTHDF2-depleted cells (Fig. 4N and O). Similarly, only wild-type YTHDF2 rescued the decreased MYC levels upon YTHDF2 depletion (Supplementary Fig. S5G and S5H). Exogenous expression of both wild-type and mutant YTHDF2 clones had minimal effect on endogenous MYC mRNA levels (Supplementary Fig. S5I). Collectively, these results establish that the effects of YTHDF2 on MYC and VEGFA transcript stability and GSC viability depend on m6A.
YTHDF2 Acts as a GSC-Specific Dependency by Preserving MYC Transcript Stability
MYC serves as a master regulator of cancer stem cell biology, including in glioblastoma (42, 44), suggesting that YTHDF2 regulates the viability and stemness of GSCs, at least in part, through MYC. Knockdown of MYC with two independent, nonoverlapping shRNAs reduced mRNA levels of VEGFA, and other genes downregulated upon YTHDF2 depletion (Fig. 5A and B). MYC overexpression rescued GSCs from YTHDF2 knockout–induced cell death (Fig. 5C and D) and partially restored mRNA levels of YTHDF2 target genes in GSCs (Fig. 5E and F). These results establish MYC as a key downstream effector of YTHDF2 in GSCs.
To determine the specificity of YTHDF2-mediated effects on MYC in GSCs, we compared the effects of YTHDF2 depletion between NSCs and GSCs using an siRNA against YTHDF2 to avoid Cas9- and viral-mediated toxicity to NSCs. YTHDF2 knockdown in NSCs did not affect MYC mRNA levels, but reduced MYC mRNA levels in GSCs (Fig. 5G–I). A time-course after actinomycin D revealed that siRNA-mediated knockdown of YTHDF2 specifically decreased MYC mRNA stability in GSCs, with little impact in NSCs (Fig. 5J). In parallel, YTHDF2 depletion reduced GSC viability without affecting NSCs (Fig. 5K). Stem cell frequency and self-renewal capacity, as assayed by ELDA, following depletion of YTHDF2 in two NSCs were not affected by YTHDF2 knockdown (Supplementary Fig. S6A and S6B). Thus, YTHDF2 represents a GSC-specific dependency that supports glioblastoma viability through GSC-specific stabilization of MYC transcripts.
IGFBP3 Is a Downstream Target of the YTHDF2–MYC Axis in GSCs
Based on the diversity of target regulation by YTHDF2, we hypothesized that molecular targets that we identified could act within a network in GSCs. IGFBP3 was one of the top downregulated genes upon YTHDF2 depletion (Fig. 3A). Prosurvival functions of IGFBP3 have been demonstrated in various cancers (43, 46), but its role in glioblastoma is not well studied. Although we detected YTHDF2 binding on IGFBP3 mRNA, YTHDF2 did not affect IGFBP3 mRNA stability (Supplementary Fig. S7A and S7B). IGFBP3 expression was downregulated following either YTHDF2 (Supplementary Fig. S3C) or MYC (Supplementary Fig. S7C) knockdown. Thus, we investigated whether IGFBP3 regulated cell viability downstream of the YTHDF2–MYC axis. Supporting the functional relevance of IGFBP3 in GSCs, IGFBP3 depletion reduced GSC viability (Fig. 6A and B) and sphere formation (Fig. 6C and D). IGFBP3 knockdown has a relatively milder effect on the viability of NSCs with a reduction in cell viability observed only on day 5 after knockdown (Supplementary Fig. S7D and S7E). IGFBP3 overexpression rescued GSCs from YTHDF2 downregulation–mediated cell death (Fig. 6E and F). Next, we investigated the contribution of IGFBP3 to in vivo glioma growth. Analysis of the TCGA glioblastoma dataset revealed that IGFBP3 mRNA expression was upregulated in patient tumors compared with normal brain (Fig. 6G). In validation studies, we measured the expression of IGFBP3 in a cohort of 20 glioblastomas and 20 nontumor brain tissues, observing elevated IGFBP3 mRNA expression in glioblastoma tissue (Fig. 6H). IGFBP3 expression portended poor prognosis in high-grade isocitrate dehydrogenase (IDH) wild-type tumors in the TCGA glioblastoma dataset (Fig. 6I). In the TCGA glioma patient data, IGFBP3 expression was associated with tumor grade in IDH1 wild-type tumors and was highly correlated with glucose metabolism–associated genes (Supplementary Fig. S7F). Consistent with these in silico results, IGFBP3 knockdown with either of two nonoverlapping shRNAs decreased expression of the glucose metabolism genes HK2, LDHA, ENO1, and MDH2 relative to a nontargeting control shRNA (Supplementary Fig. S7G). Rescue by expression of exogenous IGFBP3 compensated for YTHDF2 depletion on metabolism-associated gene expression (except MDH2; Supplementary Fig. S7H), indicating that IGFBP3 is a downstream effector of metabolic programs regulated by YTHDF2. The knockout and overexpression efficiency of YTHDF2 and FLAG-tagged IGFBP3 were validated by immunoblotting (Supplementary Fig. S7I).
To examine the functional contributions of IGFBP3 to glioma growth, we knocked down IGFBP3 expression in GSCs and implanted these cells in the brains of immunocompromised mice. IGFBP3 knockdown in vivo led to increased survival of tumor-bearing mice and reduced tumor size (Fig. 6J–M).
To assess the contribution of IGFBP3 in YTHDF2-mediated regulation of MYC mRNA levels and stability and to rule out a potential indirect feedback effect, we measured the MYC mRNA levels upon IGFBP3 perturbations. IGFBP3 knockdown showed no effect on MYC levels in GSCs (Supplementary Fig. S7J). Further, rescue by expression of exogenous IGFBP3 failed to compensate for YTHDF2 depletion–mediated MYC expression (Supplementary Fig. S7K). Collectively, these results establish IGFBP3 as a key downstream effector of the YTHDF2–MYC signaling axis in GSCs.
The YTHDF2–MYC–IGFBP3 Axis Promotes In Vivo Tumor Growth
To address the potential benefit of therapeutic targeting of YTHDF2 in vivo, mice bearing orthotopic xenografts from patient-derived GSCs were transduced with a nontargeting, control sgRNA (sgCONT) or one of two sgRNAs targeting YTHDF2. Knockout of YTHDF2 prolonged tumor latency, as measured by time to onset of neurological signs, and reduced tumor volumes compared with mice bearing GSCs transduced with a control sgRNA (Fig. 7A–D). YTHDF2 knockout thus impaired the in vivo tumor formation capacity of GSCs. To further understand the clinical significance and pan-cancer dependency on YTHDF2, we interrogated the Cancer Dependency Map (www.depmap.org). Across a panel of 558 cell lines profiled using whole-genome CRISPR/Cas9 screening, several cancer types showed a dependency on YTHDF2 (Supplementary Fig. S8A). YTHDF2 expression correlated with poor prognosis in the TCGA glioblastoma and CGGA datasets (Supplementary Fig. S8B and S8C). Given the ability of IGFBP3 to rescue the proliferative effects of YTHDF2 depletion in vitro, we investigated whether IGFBP3 could reverse the effects of YTHDF2 knockdown in vivo. IGFBP3 overexpression restored in vivo tumor formation capacity of YTHDF2-depleted GSCs (Fig. 7E–H).
We then sought to identify YTHDF2-based druggable therapeutic dependencies. We leveraged the Cancer Therapeutics Response Portal (CTRP) and Cancer Cell Line Encyclopedia datasets, which contain mRNA expression and drug screening data from 481 compounds in 860 cancer cell lines. Correlation of drug sensitivity (AUC) with YTHDF2 mRNA expression in brain cancer cell lines revealed that high YTHDF2 expression correlated with sensitivity to a CDK9 inhibitor, a SIRT1 inhibitor (salermide), and an IGF1/IGF1R inhibitor (linsitinib), among others (Fig. 7I). We prioritized linsitinib for further investigation based upon the role of the IGF pathway member IGFBP3 as a downstream effector of YTHDF2 in GSCs. Linsitinib treatment reduced the cell viability of a panel of GSCs in a concentration-dependent manner without affecting the viability of a panel of NSCs (Fig. 7J). Linsitinib treatment blocked the phosphorylation of IGF1R and downstream signaling molecules AKT, S6K, and S6RP in 387 and 4121 GSCs in a concentration-dependent manner (Supplementary Fig. S9). In in vivo studies, oral administration of linsitinib (50 mg/kg body weight) prolonged tumor latency as measured by time to onset of neurological signs and reduced tumor volumes compared with mice bearing intracranial xenografts of two patient-derived GSC models (Fig. 7K–N).
Collectively, these studies demonstrate that the m6A reader YTHDF2 and its downstream effector IGFBP3 are promising therapeutic targets in glioblastoma and that linsitinib serves as an effective agent with a broad therapeutic window.
Discussion
Although discovered in the early 1970s, the biological significance of mRNA modifications has only recently been appreciated, owing to the discoveries of the important regulator proteins and advances in transcriptome-wide m6A mapping techniques. Although many RNA-binding proteins have been implicated in cancer development, the functional importance of m6A modifiers in the glioblastoma hierarchy remains an open area of investigation. Here, we interrogated m6A localization and its regulators by mapping m6A-methylated RNAs in GSCs and their normal counterpart, NSCs. Identification of YTHDF2 as a GSC-specific m6A effector implicates the m6A pathway as a critical dependency in glioblastoma.
Oncogenic functions for m6A “writers” (34) and “erasers” (21) indicate that target-specific m6A modifications contribute to glioblastoma tumor maintenance. The m6A writer METTL3 is reported to have both oncogenic (17) and tumor-suppressive (47) functions in glioblastoma, highlighting the complexity of the functional role of m6A RNA modifications. Although RNA and DNA methylation appear distinctly regulated, lessons gained from DNA methylation studies caution against oversimplifying methylation events as solely oncogenic or tumor-suppressive. Methylation can have diverse effects in DNA regulation, and our results support a similarly nuanced function of RNA methylation. By overlapping the m6A profiling data with RNA-seq in our GSC models, we observed that the steady-state expression of m6A-tagged transcripts was higher than non-m6A transcripts. YTHDF2 has been reported to destabilize m6A-tagged transcripts by activating mRNA degradation pathways (29, 48, 49). Although our findings are largely consistent with the proposed role of YTHDF2 in mRNA destabilization, we also observe that YTHDF2 stabilizes certain important transcripts in GSCs, suggesting that the precise nature of the interaction between YTHDF2 and transcripts and its outcome might be controlled by additional unknown factors which might function in a cell type–specific manner. Alternatively, YTHDF2 may destabilize mRNAs that encode suppressors of other mRNAs, thus allowing YTHDF2 to stabilize certain transcripts through an indirect pathway. We also determined that YTHDF2 expression was elevated in GSCs as well as in glioblastoma tumors. Moreover, when we overlapped RNA-seq data from YTHDF2-knockout GSCs with m6A profiling data, m6A-tagged mRNAs tended to be downregulated in comparison with non–m6A-tagged mRNAs. These observations emphasize the complexity of the functional outcome of m6A–YTHDF2 interaction in regulating gene expression.
MYC is overexpressed in approximately 70% of malignancies (50) and informs poor prognosis in patients with glioblastoma. MYC is also a central regulator of GSCs (42). Direct targeting of MYC has been a challenge owing to its “undruggable” protein structure. Thus, several studies have aimed to identify regulators of MYC as an alternative therapeutic strategy for multiple cancer types including glioblastoma. We observed a cell type–specific effect of YTHDF2 on MYC where it selectively stabilizes MYC transcripts in GSCs, but not in NSCs, and establishes a GSC-specific dependency. These findings indicate the involvement of cell type–specific factors that may convert the general mRNA-destabilizing function of YTHDF2 to the mRNA-stabilizing function and warrant further investigation. Moreover, it will be important to interrogate whether certain elements in MYC and other mRNAs contribute to such functional alteration. Recent studies have found that YTHDF1 and YTHDF3, which are typically expressed at low levels, show minimal effects upon knockout due to compensation by the more abundant YTHDF2 (13). Nevertheless, YTHDF1 and YTHDF3 can partially compensate for the depletion of the more highly expressed YTHDF2 paralog (13). Depletion of all three YTHDF paralogs may produce even more robust effects on MYC expression and tumor growth suppression. Thus, YTHDF2, or the YTHDF family of proteins, may represent a target with a broad therapeutic index to perturb MYC signaling in glioblastoma.
We established IGFBP3 as an indirect target of YTHDF2, as YTHDF2 depletion downregulated IGFBP3 mRNA and protein levels without affecting its mRNA stability. YTHDF2 regulated IGFBP3 levels via MYC in GSCs. IGFBP3 is a member of the family of IGF-binding proteins, which promote several malignancies (51, 52). To address the translational potential of YTHDF2, we predicted available therapeutic compounds with efficacy in YTHDF2-overexpressing cells. We identified the IGF1/IGF1R inhibitor linsitinib as a selective inhibitor of GSC growth and tumor formation, which further validates our finding that IGFBP3 is a downstream effector of YTHDF2. Linsitinib has been used in several clinical trials and has been well tolerated (53, 54), although it has not yet been tested in patients with glioblastoma. Our results support further investigation into the clinical utility of linsitinib and other IGF1R inhibitors in glioblastoma.
Through a combination of in vitro and in vivo studies in GSCs, we elucidated the functions of m6A mediators in GSCs and identified YTHDF2 as a GSC-specific dependency that regulates glucose metabolism in GSCs through stabilization of MYC transcripts. These findings reveal new therapeutic opportunities for targeting previously undruggable, essential nodes in glioblastoma.
Methods
GSC Derivation
Glioblastoma tissues were obtained from excess surgical resection samples from patients at Duke University after review by neuropathology with written informed consent and in accordance with an Institutional Review Board–approved protocol (090401). All patient studies were conducted in accordance with the Declaration of Helsinki. GSC387 and GSC4121 were derived by our laboratory and transferred via a material transfer agreement from Duke University. GSC23 was transferred via a material transfer agreement from MD Anderson Cancer Center. All GSCs were cultured in 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 decrease the incidence of cell culture–based artifacts, patient-derived xenografts were produced and propagated as a renewable source of tumor cells for study. Short tandem repeat analyses were performed to authenticate the identity of each tumor model used in this study on a yearly basis. Cells were frozen and stored at -196°C (liquid nitrogen) when not being actively cultured. Mycoplasma testing was performed by qPCR cellular supernatants on a yearly basis. Cells were grown for fewer than 20 in vitro passages from xenografts.
Other Cell Models
The nonmalignant NSC models ENSA, HNP1, and NSC11 as well as normal human astrocytes were used in this study. ENSA (ENStem-A) is a human embryonic stem–derived neural progenitor cell (Millipore Sigma, Cat# SCC003). NSC11 is a human-induced pluripotent-derived neural progenitor cell (Alstem, Cat# hNSC11). HNP1 (STEMEZ HNP1) is a human neural progenitor cell (Neuromics, Cat# HN60001).
In Vivo Tumorigenesis
For YTHDF2 sgRNA and IGFBP3 shRNA and IGFBP3 overexpression experiments, intracranial xenografts were generated by implanting 5,000 human-derived GSCs into the right cerebral cortex of NSG mice (NOD.Cg-Prkdcscid Il2rgtm1Wjl/SzJ, The Jackson Laboratory, RRID:IMSR_JAX:005557) at a depth of 3.5 mm under an Institutional Animal Care and Use Committee–approved protocol of the University of California, San Diego (UCSD). Brains were harvested and fixed in 4% formaldehyde, cryopreserved in 30% sucrose, and then cryosectioned. Hematoxylin and eosin (H&E) staining was performed on sections for histologic analysis. In parallel survival experiments, mice were observed until the development of neurological signs. For in vivo drug treatment studies, intracranial xenografts were generated by implanting 5,000 patient-derived GSCs (387 and 4121) into the right cerebral cortex of NSG mice as described above. Mice recovered for 7 days were randomly assigned into drug versus treatment group by a blinded investigator. Mice were then treated daily with either vehicle (25 mmol/L tartaric acid) or 50 mg/kg linsitinib by oral gavage.
Healthy, wild-type male or female mice of NSG background, 4 to 6 weeks old, were randomly selected and used in this study for intercranial injection. Mice had not undergone prior treatment or procedures. Mice were maintained in 14-hour light/10-hour dark cycle by animal husbandry staff with no more than 5 mice per cage. Experimental animals were housed together. 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 weight loss. Brains were harvested and fixed in 4% formaldehyde for 48 hours, stored at 4°C in 70% ethanol, and sectioned. H&E staining was performed on sections for histologic analysis.
Plasmids and Cloning
Lentiviral constructs expressing nonoverlapping shRNAs directed against YTHDF2 (TRCN0000254410 and TRCN0000168751), MYC (TRCN0000039640 and TRCN0000039642), IGFBP3 (TRCN0000072508 and TRCN0000072512), or a nontargeting control shRNA (TRCN0000231489) with no targets in the human genome were obtained from Sigma-Aldrich. All shRNAs were assayed for knockdown efficiency by qPCR and were then used for all following experiments. For CRISPR/Cas9 experiments, sgRNAs were cloned into the lentiCRISPR v2 plasmid (Addgene #52961). sgRNAs used in this study were sgYTHDF1 #1: 5′-AGTTTCAAAGCCGACCTCGT-3′, sgYTHDF1 #3: 5′-TTGGCACCTGGGATAACAAG-3′, sgYTHDF2 #1: 5′-TGGGTAGCACCTCGGAACCG-3′, sgYTHDF2 #2: 5′-GTCCATTACTAGTAACATCG-3′, sgYTHDF3 #2: 5′-TGGTGTATTTAGTCAACCTG-3′, sgYTHDF3 #4: CAACCGAAACTTAAACCCAA-3′, sgMETTL3: 5′- ATTCTGTGACTATGGAACCA-3′, sgMETTL14 #2: 5′-GTCCAGTGTCTACAAAATGT-3′, sgMETTL14 #3: 5′-TTGAAGAATATCCTAAACTG-3′. Human YTHDF2 cDNA and a control vector were obtained from Vectorbuilder (YTHDF2: VB170515–1074cyz, control: VB171005–1094yhn). Lentiviral overexpression plasmids for human MYC and IGFBP3 were generated by cloning the respective open reading frames (ORF) with an N-terminal FLAG tag into the pCDH-MCS-T2A-Puro-MSCV vector (System Biosciences) using In-Fusion HD Cloning Kit (Clontech, 638920). All plasmids were verified by Sanger sequencing. 293T cells (ATCC Cat# CRL-3216, RRID:CVCL_0063) were used to generate lentiviral particles by cotransfecting the packaging vectors psPAX2 and pMD2.G using LipoD293 transfection reagent (SignaGen, SL100668) following the manufacturer's instructions.
Proliferation and Neurosphere Formation Assays
Cell proliferation experiments were conducted by plating cells of interest at a density of 3,000 cells per well in a 96-well plate with 8 replicates. CellTiter-Glo (Thermo Fisher Scientific) was used to measure cell viability. Data are presented as mean ± SD. Neurosphere formation was measured by in vitro limiting dilution assay, as previously reported (55). Briefly, decreasing numbers of cells per well (1, 5, 10, 20, 50) were plated into 96-well plates. The presence and number of neurospheres in each well were recorded 7 days after plating. Extreme limiting dilution analysis was performed using software available at http://bioinf.wehi.edu.au/software/elda, as previously described (55, 56).
Immunoblotting
Cells were collected and lysed in RIPA buffer (50 mmol/L Tris-HCl, pH 7.5; 150 mmol/L NaCl; 0.5% NP-40; and 50 mmol/L NaF with protease inhibitors) and incubated on ice for 1 hour. Lysates were centrifuged at 4°C for 20 minutes at 14,000 rpm, and supernatant was collected. The bicinchoninic acid (BCA) assay (Bio-Rad Laboratories, RRID:SCR_008426) was utilized for determination of protein concentration. Equal amounts of protein samples were mixed with SDS Laemmli loading buffer, boiled for 5 minutes, electrophoresed using NuPAGE Bis-Tris Gels, and then transferred onto Nitrocellulose membranes. TBS-T supplemented with 5% nonfat dry milk was used for blocking for a period of 1 hour followed by blotting with primary antibodies and appropriate dilution at 4°C for 16 hours. Blots were washed 5 times for 5 minutes each with TBS-T and then incubated with appropriate secondary antibodies in 5% nonfat milk in TBS-T for 1 hour. For all Western immunoblot experiments, blots were imaged using Bio-Rad Image Lab software (RRID:SCR_008426) and subsequently processed using Adobe Illustrator (RRID:SCR_010279) to create the figures.
Immunofluorescence Staining and Imaging
For immunofluorescence microscopy, cells were plated on Matrigel-coated coverslips. After 24 hours, cells were fixed twice using 4% paraformaldehyde (15 minutes each time) and processed as described previously (57). Briefly, fixed cells were washed in PBS, neutralized (10 minutes; 50 mmol/L glycine in PBS), blocked in PBS with 1 mg/mL BSA (blocking buffer) for 10 minutes, and permeabilized in blocking buffer containing 0.05% saponin. Cells were incubated with YTHDF2 antibody (Proteintech Cat# 24744-1-AP, RRID:AB_2687435) at a 1:200 dilution in blocking buffer at 4°C overnight. The next day, the cells were washed 3 times with blocking buffer and incubated with donkey anti-rabbit secondary antibody (Molecular Probes Cat# A-21206, RRID:AB_2535792). Cells were washed 3 times in blocking buffer, and coverslips were mounted using Prolong Gold Antifade (Life Technologies). Optical section Z-stacks were imaged using 60x Magnification on Leica Confocal SPE (Sanford Consortium, UCSD facility) and processed using ImageJ software (NIH, Bethesda, MD, RRID:SCR_003070).
Quantitative RT-PCR
TRIzol reagent (Sigma Aldrich) was used to isolate total cellular RNA from cell pellets according to the manufacturer's instructions. The high-capacity cDNA reverse transcription kit (Thermo Fisher scientific, catalog 4368814) was used for reverse transcription into cDNA. Quantitative real-time PCR was performed using Applied Biosystems 7900HT cycler using Radiant Green Hi-ROX qPCR kit (catalog number QS2050). qPCR primers used in this study were SOX2-fwd: 5′-TACAGCATGTCCTACTCGCAG-3′, SOX2-rev: 5′-GAGGAAGAGGTAACCACAGGG-3′, OLIG2-fwd: 5′-TGGCTTCAAGTCATCCTCGTC-3′, OLIG2-rev: 5′-ATGGCGATGTTGAGGTCGTG-3′, GAPDH-fwd: 5′-GGAGCGAGATCCCTCCAAAAT-3′, GAPDH-rev: 5′-ATGGCGATGTTGAGGTCGTG-3′, 18S-fwd: 5′-GGCCCTGTAATTGGAATGAGTC-3′, 18S-rev: 5′-CCAAGATCCAACTACGAGCTT-3′, MYC-fwd: 5′-GGCTCCTGGCAAAAGGTCA-3′, MYC-rev: 5′-CTGCGTAGTTGTGCTGATGT-3′, YTHDF2-fwd: 5′-TAGCCAACTGCGACACATTC-3′, YTHDF2-rev: 5′-CACGACCTTGACGTTCCTTT-3′, IGFBP3-fwd: 5′-AGAGCACAGATACCCAGAACT-3′, IGFBP3-rev: 5′-GGTGATTCAGTGTGTATTCCATT-3′, VEGFA-fwd: 5′-ATCTGCATGGTGATGTTGGA-3′, VEGFA-rev: 5′-GGGCAGAATCATCACGAAGT-3′, HK2-fwd: 5′-TGCCACCAGACTAAACTAGACG-3′, HK2-rev: 5′-CCCGTGCCCACAATGAGAC-3′, LDHA-fwd: 5′-TTGACCTACGTGGCTTGGAAG-3′, LDHA-rev: 5′-GGTAACGGAATCGGGCTGAAT-3′, MDH2-fwd: 5′-GCCATGATCTGCGTCATTGC-3′, MDH2-rev: 5′-CCGAAGATTTTGTTGGGGTTGT-3′, ENO1-fwd: 5′-TGGTGTCTATCGAAGATCCCTT-3′, ENO1-rev: 5′-CCTTGGCGATCCTCTTTGG-3′, RPS2-fwd: 5′-GGCCTCTCTCAAGGATGAGGT-3′, RPS2-rev: 5′-GTCCCCGATAGCAACAAATGC-3′, RPL13a-fwd: 5′-GCCATCGTGGCTAAACAGGTA-3′, RPL13a-rev: 5′-GTTGGTGTTCATCCGCTTGC-3′, DDIT-fwd: 5′-TGAGGATGAACACTTGTGTGC-3′, DDIT-rev: 5′-CCAACTGGCTAGGCATCAGC-3′, CBP-fwd: 5′-ACCGGTGTAAGGAAAGGCTG-3′, CBP-rev: 5′-TCAGGTGTTGGGAAGATGGC-3′.
Patient Database Bioinformatics
For survival analyses, TCGA data (58) were downloaded using the “TCGA2STAT” R package (59). The Kaplan–Meier survival analysis with the log-rank analysis was used to assess prognostic significance of every gene in the TCGA GBM HG-U133A microarray and glioblastoma (GBM) or GBM-LGG RNA-seq datasets (1). The processed UCSC TOIL analysis of TCGA and GTEx RNA-seq data was used to determine genes that were differentially expressed between glioblastoma specimens and normal brain specimens (60). The Cox Proportional Hazards model and log-rank analysis were used to assess prognostic significance of every gene in the TCGA GBM HG-U133A microarray dataset regardless of IDH mutation status.
RNA-seq and Data Analysis
TRIzol reagent (Sigma Aldrich) was used to isolate total cellular RNA from cell pellets according to the manufacturer's instructions. RNA was purified and subjected to RNA-seq. FASTQ sequencing reads were trimmed using Trim Galore (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), and transcript quantification was performed using Salmon (RRID:SCR_017036) in the quasi-mapping mode (61). Salmon “quant” files were converted using Tximport (RRID:SCR_016752; https://bioconductor.org/packages/release/bioc/html/tximport.html), and differential expression analysis was performed using DESeq2 (DESeq, RRID:SCR_000154; ref. 62). GSEA was performed by selecting differentially expressed genes (FDR-corrected P value < 0.05), generating a preranked list, and inputting a preranked list into the GSEA desktop application (http://software.broadinstitute.org/gsea/downloads.jsp; refs. 63, 64). Pathway enrichment bubble plots were generated using the Bader Lab Enrichment Map Application (65) and Cytoscape (RRID:SCR_003032; http://www.cytoscape.org).
N6-methyladenosine RIP Followed by Next-Generation Sequencing
Methyl RIP (meRIP) was performed as described previously (66, 67). Briefly, total RNA was extracted using TRIzol reagent (Invitrogen) and treated with RNase-free DNase I (Roche) to deplete DNA contamination. PolyA RNA was purified using a GenElute mRNA Miniprep Kit (Sigma-Aldrich) as per the manufacturer's instructions and fragmented using an RNA fragmentation kit (Ambion). Note that 200 μg of fragmented RNA was incubated with 3 μg anti-m6A (Synaptic Systems, Cat# 202 003) in RIP buffer (150 mmol/L NaCl, 10 mmol/L Tris, and 0.1% NP-40) for 2 hours at 4°C, followed by addition of washed protein A/G magnetic beads (Millipore) and incubation at 4°C for a further 2 hours. Beads were washed 6 times in RIP buffer and incubated with 50 μL immunoprecipitation buffer containing 0.5 mg/mL m6AMP (Sigma-Aldrich) to elute RNA. Immunoprecipitated RNA was extracted with phenol/chloroform, and RNA samples were sent for next-generation sequencing.
meRIP-seq Alignment and Peak Calling
FASTQ files were processed using Trim Galore v0.4.5 (RRID:SCR_011847) and cutadapt v2.1 (RRID:SCR_011841) to trim adaptors. FASTQ files were aligned using STAR v2.7 (RRID:SCR_015899) to hg19 and sorted and indexed using SAMtools v1.3 (RRID:SCR_002105). Indexed bams were processed to generate individual and differential peaksets using m6a Viewer v1.6. Peaks were detected using default settings comparing meRIP followed by sequencing (meRIP-seq) with input at a minimum 2-fold enrichment and FDR ≤ 0.05 using Benjamini–Hochberg correction. To generate consensus peak sets for each cell model, overlapping peaks were used. Differential peaks between NSCs and GSCs were called using the R package MeTDiff with an FDR cutoff of 0.05 and differential fold enrichment of 1. m6A metagene plots were generated using the “makeMetaGeneProfile.pl” function in HOMER.
RIP-qPCR Analysis
RIP for YTHDF2 was performed using 2 μg of anti-YTHDF2 (Cat# 24744–1-AP, Proteintech), as described earlier (68). Briefly, cell extracts from 10 × 106 cells were prepared as described, and 10% volume was taken out as input at this stage. Remaining cell extracts were incubated with anti-YTHDF2 antibody coated with magnetic Protein A/G beads for 12 to 16 hours for IP. After IP, RNA was isolated from Input and IP fractions using phenol/chloroform extraction. cDNA was prepared using high-capacity cDNA reverse transcription kit (Thermo Fisher Scientific, catalog 4368814). Quantitative real-time PCR was performed using Applied Biosystems 7900HT cycler using Radiant Green Hi-ROX qPCR Kit (catalog number QS2050).
Construction of YTHDF2 iCLIP Libraries
YTHDF2 iCLIP was performed as previously described (69, 70) with the following modifications. Cross-linked cell pellets were resuspended in 100 μL of cold 1% SDS + 10 mmol/L DTT + 10X Halt protease inhibitor cocktail (Thermo Fisher) solution to dissociate YTHDF2 complexes, and sonicated using a Branson Digital Sonifier Model 450 fitted with 3.125 mm tapered micro tip probe on ice at 10% amplitude (2-second ON, 10-second OFF cycle, total 20 seconds). The pellet was topped up to 1mL of 1X RIPA buffer (50 mmol/L Tris-HCl, pH 7.5; 200 mmol/L NaCl; 1% NP-40; 0.5% sodium deoxycholate; and 0.1% SDS) and sonicated again at 10% amplitude (2-second ON, 10-second OFF cycle, total 20 seconds). Protein concentration was measured using the Pierce BCA Protein Assay Kit. DNase I (New England Biolabs) and RNase I (Thermo Fisher) were used to digest DNA and RNA in the cell pellet (high RNase I dilution = 1:5, low RNase I dilution = 1:50 or 1:100 depending on cell line). Note that 4 μg of YTHDF2 antibody (Aviva Systems Biology, ARP67917) was used for IP of each 1 mg of cell lysate. Equivalent amounts of normal rabbit IgG (Cell Signaling Technology) were used for IP as a negative control. YTHDF2 antibodies or Rabbit IgGs were first bound to Pierce Protein A/G Magnetic Beads (Thermo Fisher) in RIPA buffer for 1 hour at 4°C. IP was performed in RIPA buffer at 4°C overnight. RNA bound to YTHDF2 was ligated to adaptors and purified as previously described (69). SuperScript IV and primers with unique barcodes were used for reverse transcription to tag each replicate. The medium cDNA fraction (80–100 nucleotides) was extracted by gel purification and sequenced on the Illumina HiSeq4000 (paired end, 150 bases). Due to low levels of YTHDF2-bound RNA in the cell line NSC11, YTHDF2 iCLIP libraries were not sequenced for the NSC11 cell line.
Analysis of YTHDF2 iCLIP Libraries
Low-quality bases and adaptor sequences were trimmed using FLEXBAR (71) to remove reads with PHRED (RRID:SCR_001017) quality score < 30. Reads were demultiplexed using the pyCRAC software's pyBarcodeFilter (72) and aligned to hg19 using the bwa alignment tool (ref. 73; bwa aln -n 0.06). Paired-end alignments were pooled as previously described to allow downstream analysis without duplicate sequence names (74). For calling sites with truncations, only forward reads were used. Duplicates were pooled for analysis using the CTK pipeline (75) for identification of CITS sites and CIMS deletion sites (FDR < 0.05). For motif analysis, CITS and CIMS deletion sites were pooled, and sites in genic regions were extracted. Motifs were called from DNA sequences ±10 nucleotides of these sites using HOMER (RRID:SCR_010881). To identify the positions of the YTHDF2 iCLIP sites along an mRNA, metagene analysis was performed on the pooled cross-linking–induced truncation sites (CITS) and cross-linking–induced mutation sites (CIMS) using metaplotR (76).
RNA Stability Assay
Actinomycin D (Sigma-Aldrich) at 5 μg/mL was either added to GSCs or transfected with YTHDF2 siRNA or transduced with CRISPR/Cas9 targeting YTHDF2 using two nonoverlapping sgRNAs or nontargeting control sgRNAs. After 0, 0.5, 1, 3, and 6 hours of incubation, cells were collected and RNAs were isolated for qPCR.
Drug Sensitivity Prediction
Therapeutic sensitivity and gene-expression data were accessed through the Cancer Therapeutics Response Portal (https://portals.broadinstitute.org/ctrp/; refs. 77–79). Gene signature scores were calculated for each cell line in the dataset using the single sample Gene Set Enrichment Analysis Projection (ssGSEAProjection) module on GenePattern (RRID:SCR_003201; https://cloud.genepattern.org). Gene signature score was then correlated with AUC values for drug sensitivity for each compound tested. Correlation r value was plotted, and statistical analyses were corrected for multiple test correction.
Quantification and Statistical Analysis
All statistical analyses are described in the figure legends. For TCGA glioblastoma versus normal brain RNA-seq calculations, four-way ANOVA controlling for sex, age, and ethnicity with Benjamini and Hochberg FDR method was used for statistical analysis. For survival analyses, the Cox proportional hazards and log-rank analyses were used. For qPCR analyses, the Student t test was used to assess statistical significance, when appropriate. Two-way repeated measures ANOVA was used for statistical analysis with Dunnett multiple hypothesis test correction. For proliferation assays, two-way repeated measures ANOVA was used for statistical analysis with Dunnett multiple hypothesis test correction.
Data and Software Availability
All newly generated raw sequencing data related to this study (related to YTHDF2-knockout RNA-seq, RNA m6A profiling in NSCs and GSCs, and iCLIP-seq) are available on the Gene Expression Omnibus using the accession number GSE158742. All data accessed from external sources and prior publications have been referenced in the text and corresponding figure legends. Additional data will be made available upon request.
Contact for Reagent and Resource Sharing
Further information and requests for resources and reagents should be directed to and will be fulfilled by corresponding author Jeremy N. Rich ([email protected]).
Authors' Disclosures
S.R. Jaffrey reports personal fees and nonfinancial support from Gotham Therapeutics (S.R. Jaffrey is scientific founder, advisor to, and owns equity in Gotham Therapeutics) outside the submitted work; in addition, S.R. Jaffrey has a patent for mapping m6A in the transcriptome pending. J.N. Rich reports grants from NIH during the conduct of the study. No disclosures were reported by the other authors.
Authors' Contributions
D. Dixit: Conceptualization, resources, data curation, formal analysis, validation, investigation, visualization, methodology, writing-original draft, writing-review and editing. B.C. Prager: Resources, software, formal analysis, visualization, methodology. R.C. Gimple: Resources, software, formal analysis, visualization, methodology. H.X. Poh: Resources, data curation, formal analysis, visualization, writing-review and editing. Y. Wang: Resources, data curation, methodology. Q. Wu: Supervision, methodology. Z. Qiu: Formal analysis. R.L. Kidwell: Data curation. L.J.Y. Kim: Data curation. Q. Xie: Data curation. K. Vitting-Seerup: Software. S. Bhargava: Data curation. Z. Dong: Data curation. L. Jiang: Data curation. Z. Zhu: Data curation. P. Hamerlik: Software, formal analysis, methodology. S.R. Jaffrey: Resources, software, formal analysis, writing-review and editing. J.C. Zhao: Conceptualization, resources, formal analysis, project administration, writing-review and editing. X. Wang: Conceptualization, data curation, visualization, methodology, writing-review and editing. J.N. Rich: Conceptualization, resources, formal analysis, supervision, funding acquisition, investigation, visualization, writing-original draft, project administration, writing-review and editing.
Acknowledgments
We appreciate Dr. Xiang-Dong Fu and members of the Rich laboratory for critical reading of the article and helpful discussions. Control and m6A-binding mutant YTHDF2 constructs were generously provided by Dr. Bin Shen from Nanjing Medical University, Nanjing, China. We appreciate the Histology Core at UCSD for their work on histologic experiments and analysis. Finally, we would like to thank our funding sources: NIH grants CA217066 (B.C. Prager), CA217065 (R.C. Gimple), CA203101 (L.J.Y. Kim), CA186702 (S.R. Jaffrey), GM110090, GM132292 (J.C. Zhao), and CA197718, CA238662, and NS103434 (J.N. Rich). P. Hamerlik is supported by the Danish Cancer Society/Kraeftens Bekaempelses Videnskabelige Udvalg Foundation, Knaek Cancer, and the Novo Nordisk Foundation. H.X. Poh is supported by Agency for Science, Technology and Research (A*STAR).
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.