Research ArticleCANCER

Cross-talk among writers, readers, and erasers of m6A regulates cancer growth and progression

See allHide authors and affiliations

Science Advances  03 Oct 2018:
Vol. 4, no. 10, eaar8263
DOI: 10.1126/sciadv.aar8263

Abstract

The importance of RNA methylation in biological processes is an emerging focus of investigation. We report that altering m6A levels by silencing either N6-adenosine methyltransferase METTL14 (methyltransferase-like 14) or demethylase ALKBH5 (ALKB homolog 5) inhibits cancer growth and invasion. METTL14/ALKBH5 mediate their protumorigenic function by regulating m6A levels of key epithelial-mesenchymal transition and angiogenesis-associated transcripts, including transforming growth factor–β signaling pathway genes. Using MeRIP-seq (methylated RNA immunoprecipitation sequencing) analysis and functional studies, we find that these target genes are particularly sensitive to changes in m6A modifications, as altered m6A status leads to aberrant expression of these genes, resulting in inappropriate cell cycle progression and evasion of apoptosis. Our results reveal that METTL14 and ALKBH5 determine the m6A status of target genes by controlling each other’s expression and by inhibiting m6A reader YTHDF3 (YTH N6-methyladenosine RNA binding protein 3), which blocks RNA demethylase activity. Furthermore, we show that ALKBH5/METTL14 constitute a positive feedback loop with RNA stability factor HuR to regulate the stability of target transcripts. We discover that hypoxia alters the level/activity of writers, erasers, and readers, leading to decreased m6A and consequently increased expression of target transcripts in cancer cells. This study unveils a previously undefined role for m6A in cancer and shows that the collaboration among writers-erasers-readers sets up the m6A threshold to ensure the stability of progrowth/proliferation-specific genes, and protumorigenic stimulus, such as hypoxia, perturbs that m6A threshold, leading to uncontrolled expression/activity of those genes, resulting in tumor growth, angiogenesis, and progression.

INTRODUCTION

In contrast to the well-established role of DNA methylation and histone modifications, the role of RNA methylation in cellular processes has just begun to capture the imagination of the scientific community. There are more than 100 posttranscriptional modifications known to occur in cellular RNAs (1, 2). Of those, N6-methyladenosine is the predominant modification of mRNAs in mammalian cells (1, 2). m6A methylation is catalyzed by the multicomponent RNA methyltransferase complex, RNA demethylases, and m6A readers. The core components of the RNA methyltransferase complex include methyltransferase-like 3 (METTL3), METTL14, and Wilms tumor 1 associated protein (WTAP) (3). WTAP is a spliceosome-associated protein, which was recently shown to interact with RNA methyltransferases METTL3 and METTL14 and recruit them to the site of methylation (3). Examples of RNA demethylases include ALKB homolog 5 (ALKBH5) and fat mass and obesity-associated protein (FTO), which act as erasers of m6A and N6,2′-O-dimethyladenosine (m6Am) marks on target transcripts (4). In addition to methyltransferase complex and demethylases, m6A status is governed by the activities of m6A readers such as YTH N6-methyladenosine RNA binding proteins (YTHDF) (5). Researchers have suggested that m6A methylation plays important roles in the regulation of gene expression by affecting RNA stability, mRNA degradation, and translation (6, 7). In addition, Batista et al. (8) have recently shown that m6A methylation plays a critical role during embryonic stem cell fate decisions. Furthermore, Jia et al. (9) proved that FTO regulates energy homeostasis. In addition, Fustin et al. (10) proposed that m6A regulates human circadian rhythm. A role for m6A methylation in cancer has also begun to emerge (11, 12). For example, Li et al. (13) have recently shown that FTO plays a role in acute myeloid leukemia. In addition, Zhang et al. (14, 15) have demonstrated that hypoxia-dependent ALKBH5 expression regulates breast cancer stem cell population. Furthermore, Cui et al. have recently reported that m6A RNA methylation regulates the self-renewal of glioblastoma stem cells (16). Despite these recent developments, much needs to be done to understand the mechanism by which oncogenic stimuli may regulate the m6A methylation levels of specific transcripts and whether altered m6A levels of those genes may play a role in tumorigenesis.

Here, we report that interplay among m6A writer, eraser, and reader determines the m6A levels and, consequently, the stability of several transcripts that are known to play a critical role in cell cycle, epithelial-mesenchymal transition (EMT), and angiogenesis. Examples of those genes included cyclin E1, cyclin D1, Cdk2, cyclin A2, transforming growth factor–β (TGFβ), SMAD3, vascular endothelial growth factor A (VEGFA), platelet-derived growth factor (PDGF), connective tissue growth factor (CTGF), and high-mobility group A2 (HMGA2). Our results reveal that increased expression of METTL14 and ALKBH5 led to uncontrolled activity of these genes, leading to aberrant cell cycle progression, evasion of apoptosis, and tumor progression. Our results reveal that METTL14 and ALKBH5 inhibit the expression of YTHDF3, which, in turn, blocks RNA demethylase activity and alters the m6A status of target transcripts in cancer cells. Furthermore, we demonstrated that ALKBH5/METTL14 regulate each other’s stability as well as the stability of target transcripts by activating HuR, which is known to stabilize mRNAs by binding to cis-acting adenylate-uridylate–rich elements (17). We found that hypoxia regulates the expression of ALKBH5, METTL14, and YTHDF3 and consequently alters the m6A modification and expression of target transcripts. Together, for the first time, this study shows that maintaining an optimal level of RNA methylation is critical for regulating the expression of key transcripts, and any stress/oncogenic trigger that alters the expression of m6A writers-readers-erasers and, consequently, m6A levels of those genes, results in their uncontrolled expression/activity, leading to tumor growth and progression.

RESULTS

M6A plays a critical role in cancer growth and progression

To understand the role of m6A in oncogenesis, we depleted the “writer” RNA methyltransferase METTL14 and “eraser” RNA demethylase ALKBH5 expression in cancer cells using two sets of small interfering RNAs (siRNAs) targeting different parts of the METTL14 and ALKBH5 transcripts. Both siRNAs resulted in significantly reduced levels (>90%) of METTL14 and ALKBH5 at the RNA as well as protein levels (fig. S1, A to I). METTL14 is known to physically and functionally interact with the RNA methyltransferase METTL3 and the RNA binding protein WTAP to constitute the m6A methyltransferase complex (1, 3). Furthermore, METTL14 and METTL3 are known to regulate each other’s stability at the protein level (18). Consistent with this finding, METTL14 knockdown (KD) resulted in reduced METTL3 protein expression, while there were no significant changes observed at the transcript levels (fig. S1, A, B, F, H, and I). METTL14 depletion also resulted in reduced levels of WTAP protein without affecting WTAP transcript levels (fig. S1, A, C, and F). These findings suggest that METTL14 depletion leads to impaired stability of the m6A methyltransferase complex. Consistent with these findings, METTL14 KD resulted in a 40 to 55% reduction in overall m6A levels (fig. S1J). Next, we determined the effect of METTL14 depletion on cancer cells. METTL14 KD (using either siRNA #1 or #2) led to reduced long-term viability, migration, and invasion of breast cancer (MDA-MB-231, MDA-MB-468, and BT-549), prostate cancer (DU-145), liver cancer (HepG2), and cervical cancer (HeLa) cells (Fig. 1, A and B, and figs. S1, K to P, and S2, A to D). Silencing of ALKBH5 resulting in increased overall m6A levels (fig. S2E) using either siRNA #1 or siRNA #2 also led to reduced long-term viability, migration, and invasion of cancer cells (Fig. 1, C and D, and fig. S2, F to M). To confirm our in vitro results, we performed tumor xenograft studies using siRNAs against METTL14 and ALKBH5 that showed maximum KD efficiency. Silencing of either METTL14 or ALKBH5 resulted in reduced tumor growth in mouse tumor xenograft models (Fig. 1, E and F). Next, we asked whether METTL14 or ALKBH5 silencing affected normal cells. Unlike cancer cells, depletion of METTL14 or ALKBH5 did not have any effect on the cell viability of human mammary epithelial cells (HMECs) and human embryonic kidney (HEK) 293 cells (fig. S2N). Together, these results suggested that RNA methylation might play an important role in tumorigenesis.

Fig. 1 RNA methyltransferase METTL14 and demethylase ALKBH5 promote growth and invasion of breast cancer cells.

Clonogenic assay on scrambled-siRNA– or METTL14-siRNA (METTL14 KD)–transfected (A) or ALKBH5-siRNA (ALKBH5 KD)–transfected (C) MDA-MB-231 cells. Bar graphs below show the number of colonies counted microscopically in 10 different fields. (B and D) Photomicrograph showing migrated (top) and invaded (bottom) MDA-MB-231 cells in scrambled (Scr) or METTL14 KD (B) or ALKBH5 KD (D) cells. Bar graphs show the number of migrated and invaded cells. The data shown are means ± SEM for at least three independent experiments. **P < 0.01; ***P < 0.001; ****P < 0.0001 versus control group, t test. (E and F) Photomicrographs showing representative tumor growth in nude mice injected with 2 × 106 scrambled-siRNA–transfected (control), METTL14-siRNA (METTL14 KD)–transfected (A), or ALKBH5-siRNA (ALKBH5 KD)–transfected (B) MDA-MB-231 cells mixed with Matrigel. Bar graphs show mean tumor volume for the control (n = 8), METTL14 KD (n = 8), and ALKBH5 KD (n = 8) groups at the end of the study on day 21 after implantation of the cells.

METTL14/ALKBH5 regulate key cell cycle– and angiogenesis-associated transcripts

To understand the mechanism by which METTL14 and ALKBH5 may promote cancer growth and progression, we performed RNA sequencing (RNA-seq) analyses on METTL14/ALKBH5-silenced breast cancer cells. Gene ontology analysis revealed that cell cycle progression, regulation of cell migration, EMT, and angiogenesis were some of the highly enriched biological processes that were altered in METTL14/ALKBH5 KD cells when compared with scrambled-siRNA–transfected cells (fig. S3). Consistent with this finding, cyclin D1 and cyclin E1, TGFβ1, SMAD3, matrix metallopeptidase 9 (MMP9), VEGF, platelet-derived growth factor α (PDGFA), CTGF, and high-mobility group AT-hook 2 (HMGA2) showed significantly reduced expression in METTL14 KD cancer cells (MDA-MB-231, MDA-MB-468, and BT-549) (Fig. 2A and fig. S4, A to C, E, and F). Many of the METTL14 target genes also showed significantly reduced expression in ALKBH5 KD cancer cells (Fig. 2A and fig. S4, A to C, E, and F). To further substantiate these findings, we determined target gene expression in breast cancer cells overexpressing METTL14 or ALKBH5. METTL14 or ALKBH5 overexpression led to increased expression of target genes compared to vector control (Fig. 2B and fig. S4D). Furthermore, tumor tissues from METTL14 and ALKBH5 KD xenografts showed significantly reduced levels of target genes compared to scrambled-transfected tumors (fig. S4G). To address the specificity, we determined the levels of EZH2 (Enhancer of zeste homolog 2), which did not show altered expression in our RNA-seq analysis on METTL14/ALKBH5-silenced breast cancer cells. Western blot analysis revealed no change in EZH2 levels in METTL14/ALKBH5-silenced breast cancer cells compared to scrambled-siRNA–transfected cells (fig. S4H).

Fig. 2 METTL14 and ALKBH5 support tumor growth and progression by targeting cell cycle– and TGFβ signaling–associated transcripts.

(A) Western blot analysis of target genes in scrambled-siRNA–transfected, METTL14-siRNA (METTL14 KD)–transfected, and ALKBH5-siRNA (ALKBH5 KD)–transfected MDA-MB-231, MDA-MB-468, and BT-549 cells using antibodies against the indicated proteins. Membranes were reprobed with β-actin, which served as a loading control. Gel photograph is representative of at least three independent experiments. (B) Western blot analysis of target genes in empty vector (Control) and METTL14 expression vector (METTL14 OE)–transfected and ALKBH5 expression vector (ALKBH5 OE)–transfected MDA-MB-231 cells using antibodies against the indicated proteins. Membranes were reprobed with β-actin, which served as a loading control. Gel photograph is representative of at least three independent experiments. Quantification of band intensities for (A) and (B) is shown in fig. S4 (C and D). (C) Histogram showing cell cycle distribution of scrambled-siRNA (Control), METTL14 KD (top), and ALKBH5KD (bottom) MDA-MB-231 cells. The data shown are means ± SEM of three samples for each treatment and represent three independent experiments. (D) Western blot analysis of scrambled-siRNA–transfected, METTL14-siRNA (METTL14 KD)–transfected, or ALKBH5-siRNA (ALKBH5 KD)–transfected MDA-MB-231 cells using antibodies against the indicated proteins. Vinculin served as a loading control. Gel photograph is representative of at least three independent experiments. PARP, poly(adenosine diphosphate–ribose) polymerase. (E) Histogram showing the number of annexin V+ apoptotic cells in scrambled-siRNA–, METTL14-siRNA (METTL14 KD)–, or ALKBH5-siRNA (ALKBH5KD)–transfected MDA-MB-231 cells. MDA-MB-231 cells were transfected with scrambled-siRNA or METTL14-siRNA/ALKBH5-siRNA for 48 hours before they were stained with propidium iodide (PI) and incubated with annexin V antibody and analyzed by flow cytometry. The data shown are means ± SEM of three samples for each experiment and represent three independent experiments. (F to H) Photomicrograph showing migrated MDA-MB-231 cells in scrambled-siRNA–transfected, METTL14-siRNA–transfected (F), and ALKBH5-siRNA–transfected (H) groups, treated with TGFβ1 recombinant protein. Bar graphs show the number of migrated MDA-MB-231 cells in METTL14 KD (G) and ALKBH5 KD (I) groups. *P < 0.05; ***P < 0.001; ****P < 0.0001 versus control group, t test.

The decreased expression of cell cycle genes and reduced cancer cell viability, as well as tumor growth in METTL14/ALKBH5 KD cells, prompted us to test whether m6A may regulate cancer growth by affecting cell cycle progression. Cell cycle analysis showed that cell growth was arrested in the G1-S phase in METTL14/ALKBH5-silenced cancer cells (Fig. 2C). Consistent with this finding, we observed up-regulation of the cell cycle inhibitor protein p27/Kip1 (Fig. 2D). To address whether cell cycle arrest resulted in apoptotic cell death, we determined the levels of cleaved PARP and performed annexin V staining, followed by fluorescence-activated cell sorting (FACS) analysis. METTL14/ALKBH5 depletion resulted in significantly increased cleaved PARP levels (Fig. 2D) and annexin V+ cells (Fig. 2E), while overexpression of either METTL14 or ALKBH5 blocked the chemotherapy drug doxorubicin-induced apoptosis in MDA-MB-231 breast cancer cells (fig. S4I). To further substantiate the cancer-specific effects, we determined the effects of METTL14 and ALKBH5 silencing on the apoptosis of HEK-293 cells. Annexin V staining, followed by FACS analysis, showed no significant difference in annexin V+ cells in METTL14- or ALKBH5-silenced HEK-293 cells compared to scrambled-siRNA–transfected cells (fig. S4J). These findings suggested that a change in m6A status leads to inappropriate cell cycle activity and evasion of apoptosis, two hallmarks of cancer growth and progression.

In addition to cell cycle–associated genes, TGFβ1 and other genes, including MMP9, PDGF, CTGF, and HMG2A, which are known to play a vital role in TGFβ-induced cancer metastasis and angiogenesis (19, 20), showed reduced expression in METTL14/ALKBH5 KD cells (Fig. 2 and fig. S4, A and B). These results suggest that m6A may be one of the critical regulators of cancer progression in general and TGFβ-dependent cancer progression in particular. Consistent with this assertion, TGFβ1 was found to be one of the highly altered upstream regulators in METTL14/ALKBH5-silenced cancer cells (fig. S5, A and B). Furthermore, TGFβ1 treatment rescued the inhibitory effect of METTL14/ALKBH5 KD on migration of breast cancer cells (Fig. 2, F to I). To further validate these results, we determined the levels of total and phosphorylated SMAD3 (pSMAD3) in METTL14/ALKBH5-silenced breast cancer cells. TGFβ treatment rescued the activation of SMAD3, as revealed by the increased ratio of pSMAD3 to total SMAD3 in METTL14/ALKBH5 KD cells compared to scrambled-siRNA–transfected cells (fig. S5, C and D).

METTL14/ALKBH5 constitute a positive feedback loop with HuR

To understand the mechanism by which RNA methylation may regulate target genes in cancer cells, we performed RNA stability analysis, as recent studies have suggested that m6A methylation may play an important role in regulating RNA stability and metabolism (6). METTL14-silenced breast cancer cells showed significantly reduced RNA stability of target genes when compared with scrambled-transfected cells (Fig. 3A). To gain further insight into the mechanism by which METTL14/ALKBH5 may regulate target gene stability, we tested the levels of the RNA binding protein HuR, which is reported to have increased binding to demethylated mRNA (18). Surprisingly, we found that silencing of either METTL14 or ALKBH5 using two sets of siRNAs resulted in significantly decreased protein levels of HuR in cancer cells (MDA-MB-231, MDA-MB-468, BT-549, and HeLa) (Fig. 3, B and C, and fig. S5, E and I). Consistent with this finding, overexpression of METTL14 or ALKBH5 induced the expression of HuR (fig. S5, F and G), indicating that HuR may be directly involved in regulating the stability of METTL14/ALKBH5 target genes. Furthermore, inspection of the 3′ untranslated region (UTR) of METTL14/ALKBH5 target genes revealed that m6A-binding motifs (RRACH) are in close proximity to HuR-binding sites (fig. S6). To further substantiate the role of HuR in regulating METTL14 target genes, we performed RIP using an antibody against HuR (HuR-RIP). METTL14 target genes cyclin D1, SMAD3, cyclin E1, and TGFβ1, VEGFA, HMGA2, and PDGFA were significantly enriched in immunoprecipitated samples when compared with the GAPDH (Fig. 3D). In agreement with our results, a meta-analysis of the ENCODE data set revealed that HuR binds to cyclin D1 3′UTR in cancer cell lines (fig. S6). Furthermore, CTGF is reported to be a HuR target gene (21). Moreover, computational prediction showed several potential HuR-binding sites in the 3′UTR of METTL14/ALKBH5 target genes (fig. S6). In addition to target genes, we observed enrichment of the HuR transcript in HuR-RIP samples (Fig. 3D). This is not surprising, given that HuR is known to autoregulate its own expression by interacting with its transcript (22). We also found significant enrichment of ALKBH5 and METTL14 transcripts in HuR-RIP, suggesting that HuR may regulate the mRNA stability of both m6A writers and erasers (Fig. 3D). To support these findings, we silenced HuR expression using two sets of siRNAs that target different regions of the HuR gene and evaluated expression of METTL14/ALKBH5 and their target genes. We found significantly reduced levels of METTL14, ALKBH5, and their target genes in HuR-depleted breast cancer cells (Fig. 3E and fig. S5J). In accordance with these results, overexpression of HuR in HEK-293T cells, which do not have high levels of HuR, resulted in significantly increased expression of METTL14 and ALKBH5 (fig. S5H). These results indicated that HuR regulates the expression of METTL14/ALKBH5 and their target genes by increasing their posttranscriptional stability. Next, we asked how METTL14/ALKBH5 might regulate HuR expression. Wang et al. (23) have reported that TGFβ1 and SMAD induce HuR transcriptional and translational activity. To further substantiate this, we treated cells with recombinant TGFβ1 and determined HuR levels in breast cancer cells. Introduction of TGFβ1 resulted in increased HuR expression (Fig. 3, F and G). Together, our findings indicate that HuR, METTL14/ALKBH5, and their target genes constitute a feedback loop to regulate each other’s expression in cancer cells.

Fig. 3 METTL14 and ALKBH5 constitute a positive feedback loop with HuR to regulate the stability of target genes.

(A) Quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysis showing stability of target genes in scrambled-siRNA– or METTL14-siRNA–transfected MDA-MB-231 cells treated with actinomycin (5 μg) for the indicated hours. Transcript levels in scrambled-transfected cells were normalized to 100% for each time point. The data shown are means ± SEM of three independent experiments (n = 3 biological replicates per experiment). (B and C) Western blot analysis of scrambled-siRNA–transfected, METTL14-siRNA–transfected (B), or ALKBH5-siRNA–transfected (C) MDA-MB-231, MDA-MB-468, BT-549, and HeLa cells using antibodies against the indicated proteins. The data shown are means ± SEM of three independent biological replicates. β-Actin and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) were used as loading controls. #, *, **, and $ symbols next to β-actin in (B) and (C) indicate the same loading control as in Fig. 5C. The same loading controls were used because gels were stripped and reprobed for different proteins. Relevant proteins are shown in different figures to maintain the flow of the results. Quantification of band intensities is shown in fig. S5I. (D) qRT-PCR showing enrichment of HuR, METTL14, ALKBH5, and their target genes in MDA-MB-231 cells subjected to RIP using antibody against HuR. The data shown are means ± SEM of six independent experiments. (E) Western blot analysis of the indicated proteins in two sets of scrambled-siRNA– or HuR-siRNA (KD #1 and KD #2)–transfected MDA-MD-231 cells. β-Actin and GAPDH served as loading controls. Gel photograph is representative of at least three independent experiments. Quantification of band intensities is shown in fig. S5J. (F and G) qRT-PCR (F) and Western blot (G) analysis showing HuR expression in MDA-MB-231 cells treated with or without recombinant TGFβ1 (rTGFβ1; 2 ng/ml) using HuR-specific primers and antibody. Bar graph in (G) represents band intensity quantified from all experiments using ImageJ software. The data shown are means ± SEM of three independent experiments. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001 versus control group, t test.

Collaboration among m6A writer-eraser-reader dictates the m6A status of target genes

Since m6A modification is reported to affect HuR binding and stabilization of its target genes (18), we sought to determine the m6A status of METTL14 target genes. To do so, we performed methyl RIP, followed by deep sequencing (MeRIP-seq) on control and METTL14 KD breast cancer cells. Analysis of MeRIP-seq results using our novel algorithms exomePeak and MeTDiff (24, 25) detected 15,981 and 17,312 m6A peaks from 6796 and 7194 m6A-containing transcripts in control and METTL14-silenced MDA-MB-231 cells, respectively (Fig. 4, A and B). METTL14 KD resulted in the appearance of 10,577 unique m6A peaks, while 8145 peaks disappeared (Fig. 4, A and B). Notably, METTL14 target genes were among the unique m6A peak (at one m6A site or at multiple m6A sites)–containing transcripts that appeared in the METTL14 KD breast cancer cells (Fig. 4, A and C, and fig. S7). In accordance with the previous studies (2628), GGAC was the most common m6A motif enriched in our detected peaks (fig. S7A). Furthermore, the m6A peaks were distributed in exons and especially enriched near the stop codon (fig. S7B and Fig. 4D). Next, we compared the number of enriched peaks (normalized to input) between METTL14 KD and the scrambled control cells. A total of 2483 peaks showed a significant decrease (hypo m6A peaks), while 2870 peaks showed a significant increase (hyper m6A peaks) in METTL14 KD cells compared to scrambled control cells. Integration of MeRIP-seq and RNA-seq data revealed that 422 and 330 genes with hypomethylated m6A peaks were significantly down-regulated (hypo-down) and up-regulated (hypo-up), respectively, while 428 and 455 genes with hypermethylated m6A peaks were significantly down-regulated (hyper-down) and up-regulated (hyper-up), respectively (adjusted P value <0.01; fig. S7C). The METTL14 target genes, including TGFβ1, SMAD3, cyclin E1, cyclin D1, MMP9, VEGFA, and HMGA2, showed hyper m6A with decreased expression in METTL14 KD cells compared to scrambled control cells (Figs. 2 and 4C, and fig. S7D). Furthermore, IPA of transcripts with altered m6A levels showed TGFβ1 as one of the top upstream regulators and RhoA and phosphatidylinositol 3-kinase (PI3K)–AKT signaling as some of the highly enriched biological pathways (Fig. 4, E and F). Notably, both RhoA and PI3K-AKT pathways are known to play important roles in cancer development and progression in general and to act as modulators of TGFβ-induced EMT, tumor progression, and angiogenesis, in particular (29). Next, we confirmed our MeRIP-seq results by performing real-time PCR on MeRIP samples from MDA-MB-231 and MCF-7 cells depleted for METTL14 using gene-specific primers. Consistent with our MeRIP-seq results, target genes showed hyper m6A (normalized to input) in METTL14 KD cells compared to control cells (Fig. 4G). Similarly, target transcripts in tumor tissues from METTL14 KD cells also showed hyper m6A with reduced expression compared to control xenografts (fig. S8, A to D). Similar to METTL14, target genes in ALKBH5-silenced breast cancer cells also showed hyper m6A with decreased expression when compared to control cells (fig. S8, E and F, and Fig. 2).

Fig. 4 m6A methylation analysis of control and METTL14-silenced breast cancer cells.

(A and B) MeRIP-seq analysis showing the number of peaks (A) and m6A peak–containing transcripts (B) identified in scrambled-siRNA (siCntrl)– and METTL14-siRNA (METTL14 KD)–transfected breast cancer cells. Common m6A-containing genes share at least one common peak, while unique m6A-containing genes share no peak between scrambled-siRNA and METTL14 KD breast cancer cells. (C) TGFβ1, CCNE1, and SMAD3 showing significantly enriched m6A peaks in METTL14 KD MDA-MB-231 cells compared to scrambled-siRNA. Top two tracks represent MeRIP and input for METTL14-siRNA–transfected MDA-MB-231 cells, while bottom two tracks represent MeRIP and input for scrambled-siRNA–transfected MDA-MB-231 cells. bp, base pairs; RefSeq, reference sequence. (D) Pie chart of m6A peak distribution showing proportion of total (top) and unique (bottom) peaks in different regions of genes in scrambled-siRNA and METTL14 KD cells. (E) Ingenuity Pathway Analysis (IPA) using m6A peak–containing genes shows TGFβ as one of the top upstream regulators. (F) Bar graphs show enriched canonical pathways derived from IPA using m6A-containing genes. (G) qRT-PCR showing m6A abundance (normalized to input) of target genes in MeRIP samples from MDA-MB-231 and MCF-7 cells transfected with scrambled-siRNA or METTL14-siRNA (METTL14 KD). The data shown are means ± SEM of three independent experiments (n = 3 biological replicates per experiment).

We show that METTL14, ALKBH5, and HuR constitute a positive feedback loop. Since levels of writers and erasers must be in a dynamic equilibrium to regulate m6A status, we wondered whether RNA demethylase and methyltransferase machineries may cross-talk. Silencing of ALKBH5 resulted in significantly reduced levels of METTL14, METTL3, and WTAP (Fig. 5A and fig. S8, G and H), while silencing of METTL14 reduced ALKBH5 levels in all cancer cell lines tested (Fig. 5A and fig. S8H). In accordance with this, overexpression of METTL14 or ALKBH5 resulted in increased expression of ALKBH5 and METTL14, respectively (Fig. 2B and fig. S4D). To gain more insight into the mechanism by which the m6A status of METTL14 target genes in cancer cells may be regulated, we focused on m6A reader proteins, as YTHDF2 was recently shown to block RNA demethylase activity to maintain the methylation status of heat shock stress-induced transcripts (30). Moreover, our RNA-seq analysis revealed that YTHDF3 and YTHDF2 (albeit less than YTHDF3) levels were significantly increased in METTL14 and ALKBH5 KD cells (fig. S3C). We further confirmed these results by determining YTHDF3/2 transcript and protein levels in METTL14/ALKBH5-silenced cancer cells (Fig. 5B and fig. S8, I and J). Next, we directly addressed whether YTHDF protein may block RNA demethylase activity, resulting in increased m6A levels of target genes in METTL14-silenced cells. An in vitro m6A-binding assay revealed that single-stranded TGFβ1 RNA containing m6A modification recruited YTHDF3, which, in turn, blocked RNA demethylase binding to the m6A-containing transcript (Fig. 6A). Supporting this finding, real-time PCR analysis on MeRIP samples showed that the number of m6A-containing transcripts and/or the number of m6A marks/enrichment of m6A sites in target genes were significantly decreased when both METTL14 and YTHDF3 were silenced compared to METTL14-only KD breast cancer cells (Fig. 6B). To further support our finding that YTHDF3 may block RNA demethylase activity to regulate m6A levels, we determined the distribution of YTHDF2. Using siRNA against YTHDF2, we showed that YTHDF2 protein was detected in both the cytoplasm and the nucleus (fig. S8K). These findings suggest that methyltransferase complex proteins, YTHDF proteins, and RNA demethylase machinery may collaborate to determine the methylation status of target genes in cancer cells. Next, we tested whether regulation of methylation status by writer-eraser-reader interplay may be directly associated with target gene expression. We observed increased levels of target genes in breast cancer cells depleted for YTHDF3 alone or METTL14 and YTHDF3 together compared to METTL14-only KD cells (Fig. 6, C and D, and fig. S8L).

Fig. 5 Cross-talk among writer, reader, and eraser determines m6A levels of target transcripts.

(A) Western blot analysis of scrambled-siRNA–transfected or METTL14-siRNA–transfected (top) or ALKBH5-siRNA–transfected (bottom) MDA-MB-231, MDA-MB-468, and BT-549 cells using antibodies against the indicated protein. Gel photographs represent results from three independent experiments. Note that METTL14 levels in ALKBH5 KD MDA-MB-231, MDA-MB-468, and BT-549 cells are shown in Fig. 4C. #, *, **, and $ symbols next to β-actin indicate the same loading control as in Fig. 4C. (B) Western blot analysis of scrambled-siRNA–transfected or METTL14-siRNA–transfected (top) or ALKBH5-siRNA–transfected (bottom) MDA-MB-231, MDA-MB-468, BT-549, and HeLa cells using antibodies against YTHDF3. Gel photographs represent results from three independent experiments. Quantification of band intensities for (A) and (B) is shown in fig. S8.

Fig. 6 m6A reader blocks RNA demethylase activity to regulate m6A levels and expression of target genes.

(A) m6A-contaning biotin-labeled TGFβ1 mRNA was incubated with 1 μg of FTO in the absence (−) or presence of 1 μg (+) or 2 μg (++) of YTHDF3, followed by mRNA pulldown and Western blot using antibodies against YTHDF3 or FTO. A portion of the sample collected before pulldown served as inputs. Gel photograph represents results from three independent experiments. (B) qRT-PCR showing m6A abundance (normalized to input) of target genes in MeRIP samples from MDA-MB-231 cells transfected with scrambled-siRNA, METTL14-siRNA (METTL14 KD), or METTL14 KD + YTHDF3-siRNA (YTHDF3 KD). The data shown are means ± SEM of three independent experiments. (C and D) qRT-PCR (C) and Western blot (D) analysis of scrambled-siRNA–, METTL14-siRNA (METTL14 KD)–, YTHDF3-siRNA (YTHDF3 KD)–, or METTL14-siRNA + YTHDF3-siRNA–transfected MDA-MB-231 cells using gene-specific primers and antibodies against the indicated proteins. The data shown in (C) are means ± SEM of four independent experiments. Gel photographs in (D) represent results from three independent experiments. Quantification of band intensities for (D) is shown in fig. S8K. (E) Western blot analysis of MDA-MB-231 cells exposed to normoxic and hypoxic conditions for 24 and 48 hours using antibodies against the indicated proteins. β-Actin served as a loading control. Gel photographs represent results from three independent experiments. Quantification of band intensities for (E) is shown in fig. S10A. (F) qRT-PCR showing TGFβ1 m6A abundance (normalized to input) in MeRIP samples (left) and TGFβ1 mRNA levels (right) from MDA-MB-231 cells exposed to normoxic and hypoxic conditions. The data shown are means ± SEM for two (for MeRIP) and three (for expression analysis) independent experiments. (G) qRT-PCR showing m6A abundance (normalized to input) of target gene in breast cancer patients (n = 10) and normal controls (n = 7) using gene-specific primers. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001 versus control, t test. ns, not significant. (H) Model showing mechanism by which pro-oncogenic trigger (such as hypoxia) may impair the cross-talk among m6A writer, eraser, reader, and effector proteins (such as HuR), resulting in aberrant target gene expression and, consequently, cancer growth and progression.

Next, we examined whether the cross-talk among m6A writer, eraser, and reader regulates cancer growth and progression. To address this, we performed rescue experiments by silencing YTHDF3 in METTL14- or ALKBH5-depleted cells. Remarkably, YTHDF3 silencing rescued the inhibitory effects of METTL14 and ALKBH5 KD on long-term growth and migration of breast cancer cells (fig. S9, A and B). To further investigate the role of writer and eraser axis in regulating cancer growth and progression, we overexpressed ALKBH5 or METTL14 in breast cancer cells lacking METTL14 or ALKBH5, respectively. Notably, ALKBH5 or METTL14 overexpression rescued the growth and migration of METTL14- or ALKBH5-silenced breast cancer cells, respectively (fig. S9, C and D). These results unveil a novel function for m6A writer-eraser-reader axis in regulating tumorigenesis.

Hypoxia alters expression of writer, reader, and eraser and target gene levels

Next, we wanted to address whether protumorigenic trigger/stimulus can alter the cross-talk among writer, eraser, and reader. Since hypoxia is one of the common features in most cancers, we exposed cancer cells to hypoxic conditions and determined the levels/activities of writer, reader, and eraser. Breast cancer cells exposed to hypoxic conditions showed increased expression of ALKBH5 and METTL14, while YTHDF3 level was significantly decreased compared to normoxic conditions (Fig. 6E and fig. S10A). Consistent with our results, a recent report showed induced ALKBH5 expression in breast cancer stem cells under hypoxic conditions (15). Next, we examined whether altered levels of ALKBH5/METTL14 and YTHDF3 will affect m6A level and expression of target transcript under hypoxic conditions. Real-time PCR analysis revealed that m6A levels of TGFβ1 were significantly lower, while its expression was significantly higher in breast cancer cells exposed to hypoxic conditions compared to normoxic conditions (Fig. 6F).

m6A status of target transcripts in breast cancers

To further confirm our findings and establish the physiological relevance, we performed MeRIP on tissues from breast cancer patients. qRT-PCR on immunoprecipitated samples showed the hypo m6A level of TGFβ1, cyclin D1, and SMAD3 genes in tumor tissues from breast cancer patients compared to normal matched controls (Fig. 6G). The wide variability of m6A levels in adjacent normal tissues compared to tumor tissues could be due to the presence of contaminating nonepithelial cells in normal tissues. To test the specificity of altered m6A levels in tumors from breast cancer patients, we determined m6A levels of the nuclear receptor subfamily 1, group D, member 2 (NR1D2) gene, which shows m6A-dependent expression and is known to play a role in circadian rhythm (30). Our result revealed that m6A levels of NRID2 did not change in breast cancer patients when compared to normal matched controls (Fig. 6G). Next, we analyzed whether the protumorigenic function of METTL14 and ALKBH5 may correlate with their expression levels in cancers. Meta-analysis of a large TCGA (The Cancer Genome Atlas) data set for breast, liver, and lung cancers, however, showed no significant difference in the levels of METTL14 in tumor tissues compared to normal adjacent controls (fig. S10B). As with METTL14, we observed no significant difference in ALKBH5 expression in tissues from breast, liver, and lung cancers compared with normal adjacent control tissues (fig. S10B). Since the number of control samples is smaller than that of tumor samples in the TCGA dataset, it is possible that ALKBH5 and METTL14 expression levels are higher in a subset of breast, liver, and lung cancers. It is also possible that expression of ALKBH5 and METTL14 at the protein level is higher in tumor samples compared to normal samples. A careful examination of the Human Protein Atlas dataset revealed that expression of both METTL14 and ALKBH5 in normal breast tissue is localized to patches of glandular and myoepithelial cells, while their expression in breast cancer is widespread and covers the whole tumor tissue (www.proteinatlas.org/). Next, we determined whether METTL14 and ALKBH5 transcript expression correlate with the overall survival of cancer patients. The overall survival of breast, lung, and liver cancer patients showed no significant correlation with the METTL14 or ALKBH5 expression levels (fig. S10B). In contrast to our finding, a recently published report showed lower METTL14 RNA expression in tumor tissues from cancer patients compared to normal tissue (31). This may be likely due to the limited number of samples tested for METTL14 expression in that study. We wondered whether combining METTL14, ALKBH5, and HuR together as a signature may have prognostic significance. Our analysis of the TCGA dataset revealed that, although HuR expression was significantly higher in breast cancer patients than in normal controls, HuR alone or combined expression of METTL14, ALKBH5, and HuR did not show any correlation with the overall survival of breast cancer patients (fig. S10, C and D). Together, our results indicate that m6A levels and, consequently, expression of specific tumor-associated genes, as well as ratio rather than the absolute levels of writers or erasers of RNA methylation per se, will have a greater impact on tumorigenesis (Fig. 6H).

DISCUSSION

Increasing evidence suggests that m6A RNA methylation may play vital roles in regulating cellular homeostasis (32). However, relatively little is known about the mechanism by which m6A levels of specific transcript are regulated and whether/how impaired m6A regulation of those transcripts may play a role in tumorigenesis. Our results suggest that altering m6A levels and, consequently, transcriptional stability and activity of key transcripts is critical for regulating cancer growth and progression. We provide evidence that the interplay among RNA methyltransferase, RNA demethylase, and m6A reader YTHDF dictates m6A levels of a set of transcripts, as ALKBH5 and METTL14 control each other’s expression and inhibit the expression of YTHDF3, which, in turn, blocks RNA demethylase activity and promotes degradation of target transcripts in cancer cells. We also found that METTL14/ALKBH5 activates HuR, which, in turn, regulates the stability of METTL14/ALKBH5 and target transcripts. Collectively, these studies reveal that m6A methylation of specific progrowth/proliferation transcripts is dynamically regulated, and any change in this tightly controlled process results in uncontrolled activity of those genes, leading to cancer growth and progression.

Using MeRIP-seq analysis and functional studies, we found that genes that are known to play crucial roles in cell cycle and cell proliferation were particularly sensitive to changes in m6A status. Examples of those genes included cyclin D1, Cdk2, cyclin A2, and cyclin E1. We went on to show that changes in m6A status led to aberrant expression of these genes, resulting in inappropriate cell cycle progression and evasion of apoptosis. Moreover, dysregulation of cyclin D1 and cyclin E1, which are two of the most significantly amplified genes in a wide range of cancers, contributes to loss of cell cycle control and, consequently, tumor growth. These findings support the notion that m6A may be a critical denominator that controls the activity of positive cell cycle regulators, and therefore cell growth and proliferation. In addition to genes involved in cell cycle, our studies suggest that TGFβ signaling–associated genes are particularly sensitive to changes in m6A levels. Since TGFβ is known to promote tumor growth and metastasis by inducing angiogenic factors and by facilitating EMT (33), it is likely that regulation of TGFβ signaling proteins by RNA methylation plays a critical role in tumor progression. Supporting this finding, we show that hypoxia, which is known to promote cancer progression, angiogenesis, and metastasis, altered the levels of m6A writer, eraser, and reader, resulting in decreased m6A and increased expression of TGFβ1 in cancer cells. Furthermore, we demonstrate that the expression of PDGF, which acts as a vital effector of TGFβ-mediated tumor progression (34), and HMGA2, which is used by TGFβ to mediate EMT and metastasis (35), is regulated by changes in m6A levels. In addition, we show that CTGF, which is an angiogenic mediator and a direct target of the TGFβ signaling pathway (36), is expressed at a significantly reduced level in METTL14 and ALKBH5 KD cells. m6A-dependent expression of cyclin D1 may also provide additional evidence for the importance of RNA methylation–dependent regulation of TGFβ signaling in tumor progression, as cyclin D1 is reported to play a crucial role in TGFβ-mediated breast cancer metastasis (37). In addition, our IPA of unique m6A peak–containing genes revealed that RhoA and PI3K-AKT signaling pathways, which are proven mediators of TGFβ-induced cancer progression and angiogenesis (29), were highly enriched pathways in METTL14-silenced breast cancer cells. It is worth mentioning that although several TGFβ effectors, including PDGF, CTGF, and cyclin D1, showed reduced expression in both METTL14 and ALKBH5 KD cells, HMGA2 and VEGFA levels were significantly reduced only in METTL14 KD cells. It is possible that m6A, and consequently, expression levels of these two genes, are dependent on METTL14 and some factors other than ALKBH5. Nevertheless, our results underscore that the TGFβ signaling pathway is tightly regulated by RNA methylation, and this regulation plays an important role in tumorigenesis.

Epitranscriptomic analysis showed hypo m6A for METTL14 target genes in breast cancer patients compared to normal controls. This finding is consistent with our result showing abundance of hyper m6A peaks of the target genes in METTL14 KD breast cancer cells, the epitranscriptome of which should be comparable with normal breast tissues. These results raise a pertinent question: How are m6A levels of these target transcripts regulated in cancers? Our results suggest that the ability of m6A readers to block the accessibility of RNA demethylase may be one of the critical mechanisms that determine the m6A status of target transcripts. It is plausible that the relative ratio and activity of RNA demethylases and RNA methyltransferase complex proteins determine the level of and rate at which m6A readers may occupy m6A sites on target transcripts, and in turn, regulate RNA demethylase activity. It is likely that reduced YTHDF3/2 occupancy at m6A sites, and consequently, increased RNA demethylase activity leading to decreased m6A level, may be one of the reasons for increased expression/activity of METTL14 target genes in breast cancer patients. It is tempting to speculate that the addition or removal of m6A mark(s) at certain m6A site(s) (due to combined action of writer-reader-eraser) may play a greater role in determining the stability of target transcripts. Examples of those m6A sites in target genes will include sites where unique m6A peaks appeared in METTL14 KD cells. Since m6A methylation affects mRNA stability by preventing the HuR binding to the 3′UTR of the target transcripts (18), it is also plausible that the erasure of m6A methylation at the critical m6A site(s) (due to reduced YTHDF and therefore increased RNA demethylase activities) may facilitate HuR binding and stabilization of target transcripts, leading to their increased expression/activity in cancers. Supporting this finding, we show that hypoxia, which is reported to induce HuR level/activity (38), elevates ALKBH5/METTL14 expression and inhibits YTHDF3 levels, resulting in decreased m6A and increased expression of ALKBH5/METTL14 target genes in cancer cells. Furthermore, HuR expression is shown to be elevated in breast cancer cell lines and patients when compared to normal healthy tissues (39).

Our results suggest that the m6A status of target transcripts is delicately balanced by the cross-talk/interplay among writer, eraser, and reader. It is likely that the activity/function of other RNA demethylases and other proteins with RNA methyltransferase–like activity may affect the overall number of m6A-containing transcripts and/or the number of m6A marks in m6A-containing transcripts by compensating for the loss of ALKBH5 and METTL14 in METTL14 KD or ALKBH5 KD cells, respectively. Examples of those m6A-containing transcripts will include transcripts in which m6A peaks disappeared following depletion of METTL14 KD (Fig. 4). Those m6A sites are not likely to be regulated by the YTHDF3-mediated blockage of RNA demethylase activity (as opposed to genes where unique m6A peaks appeared), either because genes containing those sites are not the target substrates for YTHDF3 or because some other factors prohibit YTHDF3 binding and blocking of RNA demethylase activity. Our results showing the importance of m6A readers in linking m6A methylation to mRNA stability and increased expression of m6A readers in METTL14/ALKBH5-depleted cells suggest that the levels of YTHDF proteins must also be dynamically regulated. It is possible that METTL14/ALKBH5 may regulate YTHDF3/2 levels by affecting the stability of co-repressors that bind to YTHDF2/3 promoters. Yet another possibility is that factors regulated by ALKBH5 and METTL14 may affect the translational stability of YTHDF3/2, which, in turn, may autoregulate their transcription. Supporting this finding, autoregulatory mechanisms are common in RNA binding proteins, as they have the ability to bind to their own mRNAs (40). Furthermore, phosphorylation of RNA binding proteins is known to affect their mRNA binding affinity (41).

In summary, adding to the increasing importance of RNA epigenetics in mammalian physiology, our study shows a critical role for m6A in tumorigenesis. Our results suggest that the cross-talk among the writer, reader, and eraser of RNA methylation may maintain a level of m6A methylation, at which binding of HuR or other RNA binding proteins to a subset of target genes, which are vital for cell proliferation, is not blocked. However, pro-oncogenic trigger/stress (for example, hypoxia) resulting in altered expression/activity of m6A machinery proteins and m6A levels may lead to unabated expression/activity of those progrowth/proliferation target genes, resulting in tumor growth and progression.

MATERIALS AND METHODS

Animals and tumor xenograft study

Athymic nude mice (6 weeks old) were obtained from Envigo Inc. Housing, and all experimental animal procedures were approved by the Institutional Animal Care and Research Advisory Committee of the University of Texas Health Science Center at San Antonio (UTHSCSA), USA. For tumor xenograft studies, MDA-MB-231 cells transfected with scrambled-siRNA or METTL14-siRNA or ALKBH5-siRNA (2 × 106) were mixed with Matrigel and injected subcutaneously in the flank of 6-week-old female athymic nude mice. Tumor volumes and body weight were measured once a week. After 21 days, mice were euthanized, and the tumors were isolated and processed for molecular studies. Tumor volume was calculated using the formula 0.5236L1 (L2)2, where L1 is the long axis and L2 is the short axis of the tumor. The number of mice used for tumor xenogaft studies are as follows: scrambled-siRNA (n = 8), METTL14-siRNA (n = 8), and ALKBH5-siRNA (n = 8).

Human cancer cells

Human cancer cell lines MDA-MB-231, MDA-MB-468, BT-549, MCF-7, HeLa, DU-145, and HepG2 were purchased from the American Type Culture Collection and cultured in standard growth medium according to their guidelines.

Migration, invasion, and colony formation assays

Cancer cells were transfected with either scrambled-siRNA or two different METTL14-siRNA/ALKBH5-siRNAs (Sigma-Aldrich) or, for 48 hours, harvested and subjected to long-term clonogenic, migration, and invasion assays, as described previously (42). For rescue experiments, MDA-MB-231 cells (2 × 105) were transiently transfected with METTL14-siRNA for 48 hours and were subsequently treated with TGFβ1 recombinant protein (2 ng/ml) for 16 hours and subjected to migration assays, as described above.

Cell cycle analysis and apoptosis assay

Both cell cycle distribution and annexin V/PI-positive cells were analyzed using flow cytometry, as described previously (42).

RNA and protein

Total RNA extracted from cancer cell lines and breast cancer and normal tissue samples were subjected to qRT-PCR and Western blot analysis. Cancer cell lines were transfected with scrambled-siRNA or METTL14-siRNA/ALKBH5-siRNA (Sigma-Aldrich) for 48 or 72 hours before they were subjected to qRT-PCR or Western blot analysis, as described previously (42). Primers used in this study are mentioned in table S1. Antibodies against ALKBH5 (#HPA 007196), β-actin (#A3854), GAPDH (#G9295), and METTL14 (#HPA038002) were purchased from Sigma-Aldrich. Antibodies against TGFβ1 (#ab155264), m6A (#ab151230), and YTHDF3 (#ab103328) were purchased from Abcam Inc. Antibodies against cyclin D1 (#2978), CDK4 (#12790), PARP (#5625S), SMAD3 (#9513), P27/KIP1 (#2552s), and pSMAD3 (#9520) were purchased from Cell Signaling Technology. The antibodies against cyclin E1 (#SC-4976), HuR (3A2) (#SC-5261), and pSMAD3 (#SC-11769) were purchased from Santa Cruz Biotechnology. The antibody against HMGA2 (GTX629478) was purchased from GeneTex. Actinomycin D (catalog no. A1410) was purchased from Sigma-Aldrich.

RNA sequencing

For whole-genome transcriptome profiling, four libraries were generated from total RNA from scrambled-siRNA– and METTL14-siRNA–transfected MDA-MB-231 cells using TruSeq RNA Library Preparation v2 according to the manufacturer’s protocol (Illumina Inc.). Samples were sequenced on the Illumina HiSeq 2000 platform (Illumina Inc.) using the 50–base pair single-end sequencing module. Sequence reads were mapped to the University of California, Santa Cruz (UCSC) genome build hg19 using TopHat aligner (v2.0.6). Gene expression was measured by HTSeq, and differential analysis was performed using DESeq. Up-regulated and down-regulated genes were determined by the following criteria: (i) absolute log2 fold change >1, (ii) average RPKM (reads per kilobase per million mapped reads) >1, and (iii) false discovery rate (FDR) <0.01. The differentially expressed genes are provided in tables S2 and S3. Raw data were deposited under Gene Expression Omnibus accession number GSE81164 (www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=ohudsiiqnjsjpmn&acc=GSE81164).

m6A MeRIP

Total RNA was isolated from scrambled and METTL14-siRNA or ALKBH5-siRNA using the RNeasy Midi Kit (Qiagen), followed by polyadenylation enrichment using the GenElute mRNA Miniprep Kit (Sigma-Aldrich). The poly (A) RNA was fragmented to 100 to 200 nucleotides using RNA fragmentation buffer (Thermo Fisher Scientific). We used 0.5 μg of fragmented mRNA as input control for RNA-seq. We used 5 μg of fragmented mRNA for m6A-containing mRNA enrichment using the Magna MeRIP m6A kit (catalog no.17-10499, Millipore), and the RNA was purified according to the manufacturer’s protocol. The resultant final product was used for RNA-seq.

MeRIP-seq and data preprocessing

The precipitated RNA was resuspended in H2O and used for library generation with the TruSeq mRNA library prep kit (Illumina). Equal amounts of barcoded samples were pooled together and submitted for cluster generation on an Illumina cBot Cluster Station and then sequencing on an Illumina HiSeq 3000 system according to the manufacturer’s instructions. The obtained sequence reads were first mapped using the TopHat aligner (v2.0.6) to the UCSC human genome build hg19. Total reads and mapped reads are listed in Table 1.

Table 1 Number of reads generated for samples by MeRIP-seq.

View this table:

Global and differential m6A methylation determination from MeRIP-seq

Differential m6A methylation analysis from MeRIP-seq data was carried out using the exomePeak software that we recently developed (24, 25). exomePeak is a freely available bioconductor package specifically designed for global m6A and differential m6A methylation analysis from MeRIP-seq data. It addresses the unique problems associated with MeRIP-seq data, including transcript abundance and transcript isoforms. To perform differential analysis, exomePeak uses a rescaled hypergeometric test. The output of exomePeak includes the loci of the differential m6A peaks, the gene symbol of the transcripts to which sites localize, the detection P values, FDRs, and the log fold change of methylation (MeFC), which is calculated for each site asEmbedded Imagewhere x and y are the read counts of the input and m6A samples, respectively, that are aligned inside a methylation site, and n and m are the total numbers of mapped reads of the respective input and m6A libraries, where a subscript 0 or 1 indicates untreated or treated (METTL14 KD) sample, respectively. MeFC directly reflects the fold change of the percentages of methylated fragments for a specific mRNA locus, where a positive value indicates high m6A methylation in the treated condition compared with the untreated condition, whereas a negative value indicates lower methylation. The output of exomePeak is then fed into the Guitar software (43), which plots the distribution of the differential m6A sites in transcripts.

In vitro m6A pulldown assay

Pulldown assay was carried out as previously described (30) with slight modification. Briefly, 100 pmol of synthesized biotinylated TGFβ mRNA containing single m6A (from Sigma-Aldrich Inc.) was incubated with streptavidin magnetic beads (Dynabeads M-280 Streptavidin, Thermo Fisher Scientific) in RNA Capture Buffer [20 mM tris (pH 7.5), 1 M NaCl, and 1 mM EDTA] at room temperature for 30 min with rotation. The bead-bound TGFβ RNA and FTO or YTHDF3 recombinant protein binding was carried out in protein RNA binding buffer [20 mM tris (pH 7.5), 50 mM NaCl, 2 mM MgCl2, and 0.1% Tween 20 detergent] at 4°C for 30 min with rotation. Beads were then washed three times with wash buffer [20 mM tris (pH 7.5), 10 mM NaCl, and 0.1% Tween 20]. 2× Laemmli sample buffer was added and boiled for 10 min to dissociate the protein from the beads. The eluted protein was detected by Western blot. A portion of the sample before pulldown was used as an input.

RNA methylation quantification

Total RNA and poly (A)+ RNA methylation were quantified using the m6A RNA Methylation Colorimetric Quantification Kit from Abcam (#ab185912) and EpiGentek (EpiQuik m6A RNA Methylation Quantification Kit; #P-9005), respectively, as per the manufacturer’s protocol. Relative m6A RNA methylation status was then calculated using the manufacturer’s supplied formula. Poly (A)+ RNA was isolated using the Dynabeads mRNA Direct Purification Kit (#61012, Life Technologies).

Statistical analysis

All values and error bars in graphs are means ± SEM. Respective n values are indicated in the figure legends. P values are determined by two-tailed Student’s t tests (GraphPad Software Inc.) and analysis of variance (ANOVA).

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/10/eaar8263/DC1

Table S1. Primer pairs and sequences used in the study.

Table S2A. Differentially expressed genes derived from comparing METTL14 KD cells versus scrambled-siRNA (control).

Table S2B. Differentially expressed genes that are highly expressed derived from comparing METTL14 KD cells versus scrambled-siRNA (control).

Table S3. Differentially expressed genes derived from comparing ALKBH5 KD cells versus scrambled-siRNA (control) using Illumina whole-genome gene expression microarrays [three microarrays for each condition (total of six microarrays) were performed for this comparison].

Table S4A. Gene set enrichment using DAVID (https://david.ncifcrf.gov/; v6.7, 2015) with 744 differentially expressed genes obtained from METTL14 KD gene expression profiling.

Table S4B. Gene set enrichment using DAVID with 440 differentially expressed genes obtained from ALKBH5 KD gene expression profiling.

Table S5A. Upstream regulators predicted by the Ingenuity Pathway Analysis (www.ingenuity.com) software with 744 DEGs of METTL14 KD gene expression profiling.

Table S5B. Upstream regulators predicted by the Ingenuity Pathway Analysis (www.ingenuity.com) software with 440 differentially expressed genes of ALKBH5 KD gene expression profiling.

Fig. S1. Efficient KD of methyltransferase complex proteins and ALKBH5 inhibits cell viability and invasion of cancer cells.

Fig. S2. METTL14 and ALKBH5 promote growth and progression of cancer cells without affecting the viability of normal cells.

Fig. S3. Cancer-associated genes are differentially expressed in METTL14/ALKBH5-silenced breast cancer cells.

Fig. S4. METTL14 and ALKBH5 regulate expression of genes involved in cell cycle, EMT, and angiogenesis.

Fig. S5. METTL14 and ALKBH5 regulate TGFβ1 and HuR expression.

Fig. S6. HuR-binding sites and m6A motif (RRACH) in 3′UTRs of METTL14/ALKBH5 target genes.

Fig. S7. Transcriptome-wide MeRIP-seq analysis shows m6A peaks in target transcripts.

Fig. S8. METTL14 and ALKBH5 regulate m6A levels of target genes by constituting a positive feedback loop and inhibiting YTHDF3.

Fig. S9. ALKBH5-YTHDF3 and METTL14-YTHDF3 axes regulate growth and migration of cancer cells.

Fig. S10. METTL14 and ALKBH5 do not show significantly different expression and association with overall survival in cancer patients.

References (44, 45)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank UTHSCSA Genome Sequencing and Greehey Children’s Cancer Research Institute Bioinformatics core facilities for performing RNA-seq and bioinformatics analysis, respectively. We thank H.-I H. Chen for TCGA data analysis. We thank L. Penalva at UTHSCSA for HuR antibody. We thank Accixx Biomedical Consulting for editing. Funding: M.K.R. was supported by NIH [National Cancer Institute (NCI)] grant R01CA179120-01A1. Y.C. was supported by NIH Clinical and Translational Science Award UL1TR001120-05 and Cancer Prevention Research Institute of Texas (CPRIT) Core Grant RP16073. Genome Sequencing Facility was supported by NIH S10OD021805-01 (to Z.L.). Y.H., Y.C., and M.K.R. were supported by NIGMS (National Institute of General Medical Sciences) R01GM113245-01 and NCI U54CA113001-10 grants. Y.C., M.K.R., and T.H.-M.H. were supported by NCI P30 grant CA054174. B.C.O. was supported by CPRIT training grant RP140105. Y.K.G. was supported by the University of Texas Rising STARs Award. P.Y. was supported by CPRIT training grant RP170345. Author contributions: S.P., V.K.E., Y.C., Y.H., and M.K.R. designed the experiments. S.P., V.K.E., P.Y., S.T., S.R., S.V., N.A., B.C.O., X.C., Z.L., and T.A.M. performed the experiments. S.P., V.K.E., Y.C., Y.H., Y.K.G., T.H.-M.H., and M.K.R. analyzed the data. S.P., V.K.E., Y.C., Y.H., and M.K.R. wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article